Purpose

Explore and analyze the continuous chlorophyll data to be included in the NDFS synthesis manuscript. We will attempt to fit multiple models to the data set and use a model selection process to determine the best one.

Global code and functions

# Load packages
library(tidyverse)
library(scales)
library(knitr)
library(mgcv)
library(car)
library(gratia)
library(here)
library(conflicted)

# Source functions
source(here("manuscript_synthesis/src/global_functions.R"))

# Declare package conflict preferences 
conflicts_prefer(dplyr::filter())

Display current versions of R and packages used for this analysis:

devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.2.3 (2023-03-15 ucrt)
##  os       Windows 10 x64 (build 19044)
##  system   x86_64, mingw32
##  ui       RTerm
##  language (EN)
##  collate  English_United States.utf8
##  ctype    English_United States.utf8
##  tz       America/Los_Angeles
##  date     2023-08-31
##  pandoc   3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version date (UTC) lib source
##  abind         1.4-5   2016-07-21 [1] CRAN (R 4.2.0)
##  bslib         0.4.2   2022-12-16 [1] CRAN (R 4.2.2)
##  cachem        1.0.8   2023-05-01 [1] CRAN (R 4.2.3)
##  callr         3.7.3   2022-11-02 [1] CRAN (R 4.2.2)
##  car         * 3.1-2   2023-03-30 [1] CRAN (R 4.2.3)
##  carData     * 3.0-5   2022-01-06 [1] CRAN (R 4.2.1)
##  cli           3.6.1   2023-03-23 [1] CRAN (R 4.2.3)
##  colorspace    2.1-0   2023-01-23 [1] CRAN (R 4.2.2)
##  conflicted  * 1.2.0   2023-02-01 [1] CRAN (R 4.2.2)
##  crayon        1.5.2   2022-09-29 [1] CRAN (R 4.2.1)
##  devtools      2.4.5   2022-10-11 [1] CRAN (R 4.2.1)
##  digest        0.6.31  2022-12-11 [1] CRAN (R 4.2.2)
##  dplyr       * 1.1.2   2023-04-20 [1] CRAN (R 4.2.3)
##  ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.2.1)
##  evaluate      0.21    2023-05-05 [1] CRAN (R 4.2.3)
##  fansi         1.0.4   2023-01-22 [1] CRAN (R 4.2.2)
##  fastmap       1.1.1   2023-02-24 [1] CRAN (R 4.2.2)
##  forcats     * 1.0.0   2023-01-29 [1] CRAN (R 4.2.2)
##  fs          * 1.6.2   2023-04-25 [1] CRAN (R 4.2.3)
##  generics      0.1.3   2022-07-05 [1] CRAN (R 4.2.1)
##  ggplot2     * 3.4.2   2023-04-03 [1] CRAN (R 4.2.3)
##  glue          1.6.2   2022-02-24 [1] CRAN (R 4.2.1)
##  gratia      * 0.8.1   2023-02-02 [1] CRAN (R 4.2.2)
##  gtable        0.3.3   2023-03-21 [1] CRAN (R 4.2.3)
##  here        * 1.0.1   2020-12-13 [1] CRAN (R 4.2.1)
##  hms           1.1.3   2023-03-21 [1] CRAN (R 4.2.3)
##  htmltools     0.5.5   2023-03-23 [1] CRAN (R 4.2.3)
##  htmlwidgets   1.6.2   2023-03-17 [1] CRAN (R 4.2.3)
##  httpuv        1.6.9   2023-02-14 [1] CRAN (R 4.2.2)
##  jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.2.1)
##  jsonlite      1.8.4   2022-12-06 [1] CRAN (R 4.2.2)
##  knitr       * 1.42    2023-01-25 [1] CRAN (R 4.2.2)
##  later         1.3.0   2021-08-18 [1] CRAN (R 4.2.1)
##  lattice       0.21-8  2023-04-05 [1] CRAN (R 4.2.3)
##  lifecycle     1.0.3   2022-10-07 [1] CRAN (R 4.2.1)
##  lubridate   * 1.9.2   2023-02-10 [1] CRAN (R 4.2.2)
##  magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.2.1)
##  Matrix        1.5-4   2023-04-04 [1] CRAN (R 4.2.3)
##  memoise       2.0.1   2021-11-26 [1] CRAN (R 4.2.1)
##  mgcv        * 1.8-42  2023-03-02 [1] CRAN (R 4.2.2)
##  mime          0.12    2021-09-28 [1] CRAN (R 4.2.0)
##  miniUI        0.1.1.1 2018-05-18 [1] CRAN (R 4.2.1)
##  munsell       0.5.0   2018-06-12 [1] CRAN (R 4.2.1)
##  mvnfast       0.2.8   2023-02-23 [1] CRAN (R 4.2.2)
##  nlme        * 3.1-162 2023-01-31 [2] CRAN (R 4.2.3)
##  patchwork     1.1.2   2022-08-19 [1] CRAN (R 4.2.1)
##  pillar        1.9.0   2023-03-22 [1] CRAN (R 4.2.3)
##  pkgbuild      1.4.0   2022-11-27 [1] CRAN (R 4.2.2)
##  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.2.1)
##  pkgload       1.3.2   2022-11-16 [1] CRAN (R 4.2.2)
##  prettyunits   1.1.1   2020-01-24 [1] CRAN (R 4.2.1)
##  processx      3.8.1   2023-04-18 [1] CRAN (R 4.2.3)
##  profvis       0.3.7   2020-11-02 [1] CRAN (R 4.2.1)
##  promises      1.2.0.1 2021-02-11 [1] CRAN (R 4.2.1)
##  ps            1.7.5   2023-04-18 [1] CRAN (R 4.2.3)
##  purrr       * 1.0.1   2023-01-10 [1] CRAN (R 4.2.2)
##  R6            2.5.1   2021-08-19 [1] CRAN (R 4.2.1)
##  Rcpp          1.0.10  2023-01-22 [1] CRAN (R 4.2.2)
##  readr       * 2.1.4   2023-02-10 [1] CRAN (R 4.2.2)
##  remotes       2.4.2   2021-11-30 [1] CRAN (R 4.2.1)
##  rlang         1.1.1   2023-04-28 [1] CRAN (R 4.2.3)
##  rmarkdown     2.21    2023-03-26 [1] CRAN (R 4.2.3)
##  rprojroot     2.0.3   2022-04-02 [1] CRAN (R 4.2.1)
##  rstudioapi    0.14    2022-08-22 [1] CRAN (R 4.2.1)
##  sass          0.4.6   2023-05-03 [1] CRAN (R 4.2.3)
##  scales      * 1.2.1   2022-08-20 [1] CRAN (R 4.2.1)
##  sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.2.1)
##  shiny         1.7.4   2022-12-15 [1] CRAN (R 4.2.2)
##  stringi       1.7.12  2023-01-11 [1] CRAN (R 4.2.2)
##  stringr     * 1.5.0   2022-12-02 [1] CRAN (R 4.2.2)
##  tibble      * 3.2.1   2023-03-20 [1] CRAN (R 4.2.3)
##  tidyr       * 1.3.0   2023-01-24 [1] CRAN (R 4.2.2)
##  tidyselect    1.2.0   2022-10-10 [1] CRAN (R 4.2.1)
##  tidyverse   * 2.0.0   2023-02-22 [1] CRAN (R 4.2.2)
##  timechange    0.2.0   2023-01-11 [1] CRAN (R 4.2.2)
##  tzdb          0.4.0   2023-05-12 [1] CRAN (R 4.2.3)
##  urlchecker    1.0.1   2021-11-30 [1] CRAN (R 4.2.1)
##  usethis       2.1.6   2022-05-25 [1] CRAN (R 4.2.1)
##  utf8          1.2.3   2023-01-31 [1] CRAN (R 4.2.2)
##  vctrs         0.6.2   2023-04-19 [1] CRAN (R 4.2.3)
##  withr         2.5.0   2022-03-03 [1] CRAN (R 4.2.1)
##  xfun          0.39    2023-04-20 [1] CRAN (R 4.2.3)
##  xtable        1.8-4   2019-04-21 [1] CRAN (R 4.2.1)
##  yaml          2.3.7   2023-01-23 [1] CRAN (R 4.2.2)
## 
##  [1] C:/R/win-library/4.2
##  [2] C:/Program Files/R/R-4.2.3/library
## 
## ──────────────────────────────────────────────────────────────────────────────

Import Data

# Define file path for processed data
fp_data <- here("manuscript_synthesis/data/processed")

# Import daily average chlorophyll and water temperature data
df_chla_wt <- readRDS(file.path(fp_data, "chla_wt_daily_avg_2013-2019.rds"))

# Import daily average flow data
df_flow <- readRDS(file.path(fp_data, "flow_daily_avg_2013-2019.rds"))

Prepare Data

# Create a vector for the factor order of StationCode
sta_order <- c(
  "RCS",
  "RD22",
  "I80",
  "LIS",
  "STTD",
  "LIB",
  "RVB"
)

# We will use LIS flow data as a proxy for STTD and RD22 flow data as a proxy
  # for I-80
df_flow_c <- df_flow %>% 
  filter(StationCode %in% c("RD22", "LIS")) %>% 
  mutate(StationCode = case_match(StationCode, "RD22" ~ "I80", "LIS" ~ "STTD")) %>% 
  bind_rows(df_flow)

# Prepare chlorophyll and flow data for exploration and analysis
df_chla_wt_c1 <- df_chla_wt %>% 
  # Join flow data to chlorophyll data
  left_join(df_flow_c, by = join_by(StationCode, Date)) %>% 
  mutate(
    # Scale and log transform chlorophyll values
    Chla_log = log(Chla * 1000),
    # Add Year and DOY variables
    Year = year(Date),
    DOY = yday(Date)
  ) %>% 
  # Add Flow Action Periods
  ndfa_action_periods() %>% 
  mutate(
    # Apply factor orders to FlowActionPeriod and StationCode
    FlowActionPeriod = factor(FlowActionPeriod, levels = c("Before", "During", "After")),
    StationCode = factor(StationCode, levels = sta_order),
    # Add a column for Year as a factor for the model
    Year_fct = factor(Year)
  ) %>% 
  arrange(StationCode, Date)

Explore sample counts by Station

df_chla_wt_c1 %>% 
  summarize(
    min_date = min(Date),
    max_date = max(Date),
    num_samples = n(),
    .by = c(StationCode, Year, FlowActionPeriod)
  ) %>% 
  arrange(StationCode, Year, FlowActionPeriod) %>% 
  kable()
StationCode Year FlowActionPeriod min_date max_date num_samples
RCS 2014 Before 2014-07-25 2014-09-08 46
RCS 2014 During 2014-09-09 2014-09-23 11
RCS 2014 After 2014-09-24 2014-09-24 1
RCS 2015 Before 2015-08-17 2015-08-20 4
RCS 2015 During 2015-08-21 2015-10-01 42
RCS 2015 After 2015-10-02 2015-11-10 40
RCS 2016 Before 2016-06-23 2016-07-13 20
RCS 2016 During 2016-07-14 2016-08-01 19
RCS 2016 After 2016-08-02 2016-09-16 46
RCS 2017 Before 2017-07-14 2017-08-28 46
RCS 2017 During 2017-08-29 2017-09-18 21
RCS 2017 After 2017-09-19 2017-11-03 46
RCS 2018 Before 2018-07-13 2018-08-27 46
RCS 2018 During 2018-08-28 2018-09-26 30
RCS 2018 After 2018-09-27 2018-11-11 46
RCS 2019 Before 2019-07-11 2019-08-25 46
RCS 2019 During 2019-08-26 2019-09-21 27
RCS 2019 After 2019-09-22 2019-11-06 46
RD22 2013 Before 2013-08-15 2013-08-21 7
RD22 2013 During 2013-08-22 2013-10-02 42
RD22 2013 After 2013-10-03 2013-11-04 33
RD22 2014 Before 2014-07-25 2014-09-08 46
RD22 2014 During 2014-09-09 2014-09-23 15
RD22 2014 After 2014-09-24 2014-11-08 46
RD22 2015 Before 2015-07-24 2015-08-20 28
RD22 2015 During 2015-08-21 2015-10-01 42
RD22 2015 After 2015-10-02 2015-11-10 40
RD22 2016 Before 2016-06-23 2016-07-13 21
RD22 2016 During 2016-07-14 2016-08-01 19
RD22 2016 After 2016-08-02 2016-09-16 46
RD22 2017 Before 2017-07-14 2017-08-28 46
RD22 2017 During 2017-08-29 2017-09-18 21
RD22 2017 After 2017-09-19 2017-11-03 46
RD22 2018 Before 2018-07-13 2018-08-27 46
RD22 2018 During 2018-08-28 2018-09-26 30
RD22 2018 After 2018-09-27 2018-11-11 46
RD22 2019 Before 2019-07-11 2019-08-25 46
RD22 2019 During 2019-08-26 2019-09-21 27
RD22 2019 After 2019-09-22 2019-11-06 46
I80 2013 Before 2013-08-15 2013-08-21 7
I80 2013 During 2013-08-22 2013-10-02 42
I80 2013 After 2013-10-03 2013-11-01 30
I80 2014 Before 2014-07-25 2014-09-08 43
I80 2014 During 2014-09-09 2014-09-23 15
I80 2014 After 2014-09-24 2014-11-08 46
I80 2015 Before 2015-07-24 2015-08-07 15
I80 2015 During 2015-08-27 2015-10-01 36
I80 2015 After 2015-10-02 2015-11-10 40
I80 2016 Before 2016-06-23 2016-07-13 21
I80 2016 During 2016-07-14 2016-08-01 19
I80 2016 After 2016-08-02 2016-09-16 46
I80 2017 Before 2017-07-14 2017-08-28 46
I80 2017 During 2017-08-29 2017-09-18 21
I80 2017 After 2017-09-19 2017-11-03 46
I80 2018 Before 2018-07-13 2018-08-27 46
I80 2018 During 2018-08-28 2018-09-26 30
I80 2018 After 2018-09-27 2018-11-11 46
I80 2019 Before 2019-07-11 2019-08-25 46
I80 2019 During 2019-08-26 2019-09-21 27
I80 2019 After 2019-09-22 2019-11-06 46
LIS 2013 Before 2013-07-07 2013-08-21 39
LIS 2013 During 2013-08-22 2013-10-02 38
LIS 2013 After 2013-10-03 2013-11-17 34
LIS 2014 Before 2014-07-25 2014-09-08 46
LIS 2014 During 2014-09-09 2014-09-23 14
LIS 2014 After 2014-09-24 2014-11-08 46
LIS 2015 Before 2015-07-06 2015-08-20 46
LIS 2015 During 2015-08-21 2015-10-01 42
LIS 2015 After 2015-10-02 2015-11-16 46
LIS 2016 Before 2016-05-29 2016-07-13 46
LIS 2016 During 2016-07-14 2016-08-01 19
LIS 2016 After 2016-08-02 2016-09-16 46
LIS 2017 Before 2017-07-14 2017-08-28 46
LIS 2017 During 2017-08-29 2017-09-09 12
LIS 2017 After 2017-09-19 2017-11-03 46
LIS 2018 Before 2018-07-13 2018-08-27 46
LIS 2018 During 2018-08-28 2018-09-26 30
LIS 2018 After 2018-09-27 2018-11-11 46
LIS 2019 Before 2019-07-11 2019-08-25 46
LIS 2019 During 2019-08-26 2019-09-21 27
LIS 2019 After 2019-09-22 2019-11-06 46
STTD 2013 Before 2013-08-15 2013-08-21 7
STTD 2013 During 2013-08-22 2013-10-02 42
STTD 2013 After 2013-10-03 2013-11-04 33
STTD 2014 Before 2014-07-25 2014-09-08 46
STTD 2014 During 2014-09-09 2014-09-23 15
STTD 2014 After 2014-09-24 2014-11-08 29
STTD 2015 Before 2015-07-27 2015-08-20 25
STTD 2015 During 2015-08-21 2015-10-01 42
STTD 2015 After 2015-10-02 2015-11-16 46
STTD 2016 Before 2016-06-23 2016-07-13 21
STTD 2016 During 2016-07-14 2016-08-01 19
STTD 2016 After 2016-08-02 2016-09-16 46
STTD 2017 Before 2017-07-14 2017-08-28 46
STTD 2017 During 2017-08-29 2017-09-01 4
STTD 2017 After 2017-09-20 2017-09-25 6
STTD 2018 Before 2018-07-20 2018-08-27 39
STTD 2018 During 2018-08-28 2018-09-26 30
STTD 2018 After 2018-09-27 2018-10-15 19
STTD 2019 Before 2019-07-26 2019-08-25 31
STTD 2019 During 2019-08-26 2019-09-21 27
STTD 2019 After 2019-09-22 2019-11-06 46
LIB 2013 Before 2013-07-16 2013-08-21 37
LIB 2013 During 2013-08-22 2013-10-02 42
LIB 2013 After 2013-10-03 2013-11-16 45
LIB 2014 Before 2014-07-25 2014-09-08 46
LIB 2014 During 2014-09-09 2014-09-23 15
LIB 2014 After 2014-09-24 2014-11-05 43
LIB 2015 Before 2015-07-06 2015-08-20 46
LIB 2015 During 2015-08-21 2015-10-01 38
LIB 2015 After 2015-10-02 2015-11-16 46
LIB 2016 Before 2016-05-29 2016-07-13 46
LIB 2016 During 2016-07-14 2016-08-01 19
LIB 2016 After 2016-08-02 2016-09-16 46
LIB 2017 Before 2017-07-14 2017-08-28 46
LIB 2017 During 2017-08-29 2017-09-18 21
LIB 2017 After 2017-09-19 2017-11-03 46
LIB 2018 Before 2018-08-14 2018-08-27 14
LIB 2018 During 2018-08-28 2018-09-26 30
LIB 2018 After 2018-09-27 2018-11-11 46
LIB 2019 Before 2019-07-11 2019-07-30 20
LIB 2019 During 2019-09-05 2019-09-21 17
LIB 2019 After 2019-09-22 2019-11-06 39
RVB 2013 Before 2013-07-07 2013-08-21 46
RVB 2013 During 2013-08-22 2013-10-02 42
RVB 2013 After 2013-10-03 2013-11-17 46
RVB 2014 Before 2014-07-25 2014-09-07 45
RVB 2014 During 2014-09-10 2014-09-23 14
RVB 2014 After 2014-09-24 2014-11-08 46
RVB 2015 Before 2015-07-06 2015-08-20 46
RVB 2015 During 2015-08-21 2015-10-01 42
RVB 2015 After 2015-10-02 2015-11-16 46
RVB 2016 Before 2016-05-29 2016-07-13 39
RVB 2016 During 2016-07-14 2016-08-01 19
RVB 2016 After 2016-08-02 2016-09-16 37
RVB 2017 Before 2017-07-14 2017-08-28 46
RVB 2017 During 2017-08-29 2017-09-18 21
RVB 2017 After 2017-09-19 2017-11-03 46
RVB 2018 Before 2018-07-13 2018-08-27 46
RVB 2018 During 2018-08-28 2018-09-26 30
RVB 2018 After 2018-09-27 2018-11-11 46
RVB 2019 Before 2019-07-11 2019-08-25 46
RVB 2019 During 2019-08-26 2019-09-21 27
RVB 2019 After 2019-09-22 2019-11-06 46

Looking at the sample counts and date ranges, we’ll only include years 2015-2019 for the analysis. However, there are still a few Station-Year-FlowPeriod combinations with low sample counts within 2015-2019.

df_chla_wt_c2 <- df_chla_wt_c1 %>% 
  filter(Year %in% 2015:2019) %>% 
  mutate(Year_fct = fct_drop(Year_fct))

df_chla_wt_c2 %>% 
  summarize(
    min_date = min(Date),
    max_date = max(Date),
    num_samples = n(),
    .by = c(StationCode, Year, FlowActionPeriod)
  ) %>% 
  arrange(StationCode, Year, FlowActionPeriod) %>% 
  kable()
StationCode Year FlowActionPeriod min_date max_date num_samples
RCS 2015 Before 2015-08-17 2015-08-20 4
RCS 2015 During 2015-08-21 2015-10-01 42
RCS 2015 After 2015-10-02 2015-11-10 40
RCS 2016 Before 2016-06-23 2016-07-13 20
RCS 2016 During 2016-07-14 2016-08-01 19
RCS 2016 After 2016-08-02 2016-09-16 46
RCS 2017 Before 2017-07-14 2017-08-28 46
RCS 2017 During 2017-08-29 2017-09-18 21
RCS 2017 After 2017-09-19 2017-11-03 46
RCS 2018 Before 2018-07-13 2018-08-27 46
RCS 2018 During 2018-08-28 2018-09-26 30
RCS 2018 After 2018-09-27 2018-11-11 46
RCS 2019 Before 2019-07-11 2019-08-25 46
RCS 2019 During 2019-08-26 2019-09-21 27
RCS 2019 After 2019-09-22 2019-11-06 46
RD22 2015 Before 2015-07-24 2015-08-20 28
RD22 2015 During 2015-08-21 2015-10-01 42
RD22 2015 After 2015-10-02 2015-11-10 40
RD22 2016 Before 2016-06-23 2016-07-13 21
RD22 2016 During 2016-07-14 2016-08-01 19
RD22 2016 After 2016-08-02 2016-09-16 46
RD22 2017 Before 2017-07-14 2017-08-28 46
RD22 2017 During 2017-08-29 2017-09-18 21
RD22 2017 After 2017-09-19 2017-11-03 46
RD22 2018 Before 2018-07-13 2018-08-27 46
RD22 2018 During 2018-08-28 2018-09-26 30
RD22 2018 After 2018-09-27 2018-11-11 46
RD22 2019 Before 2019-07-11 2019-08-25 46
RD22 2019 During 2019-08-26 2019-09-21 27
RD22 2019 After 2019-09-22 2019-11-06 46
I80 2015 Before 2015-07-24 2015-08-07 15
I80 2015 During 2015-08-27 2015-10-01 36
I80 2015 After 2015-10-02 2015-11-10 40
I80 2016 Before 2016-06-23 2016-07-13 21
I80 2016 During 2016-07-14 2016-08-01 19
I80 2016 After 2016-08-02 2016-09-16 46
I80 2017 Before 2017-07-14 2017-08-28 46
I80 2017 During 2017-08-29 2017-09-18 21
I80 2017 After 2017-09-19 2017-11-03 46
I80 2018 Before 2018-07-13 2018-08-27 46
I80 2018 During 2018-08-28 2018-09-26 30
I80 2018 After 2018-09-27 2018-11-11 46
I80 2019 Before 2019-07-11 2019-08-25 46
I80 2019 During 2019-08-26 2019-09-21 27
I80 2019 After 2019-09-22 2019-11-06 46
LIS 2015 Before 2015-07-06 2015-08-20 46
LIS 2015 During 2015-08-21 2015-10-01 42
LIS 2015 After 2015-10-02 2015-11-16 46
LIS 2016 Before 2016-05-29 2016-07-13 46
LIS 2016 During 2016-07-14 2016-08-01 19
LIS 2016 After 2016-08-02 2016-09-16 46
LIS 2017 Before 2017-07-14 2017-08-28 46
LIS 2017 During 2017-08-29 2017-09-09 12
LIS 2017 After 2017-09-19 2017-11-03 46
LIS 2018 Before 2018-07-13 2018-08-27 46
LIS 2018 During 2018-08-28 2018-09-26 30
LIS 2018 After 2018-09-27 2018-11-11 46
LIS 2019 Before 2019-07-11 2019-08-25 46
LIS 2019 During 2019-08-26 2019-09-21 27
LIS 2019 After 2019-09-22 2019-11-06 46
STTD 2015 Before 2015-07-27 2015-08-20 25
STTD 2015 During 2015-08-21 2015-10-01 42
STTD 2015 After 2015-10-02 2015-11-16 46
STTD 2016 Before 2016-06-23 2016-07-13 21
STTD 2016 During 2016-07-14 2016-08-01 19
STTD 2016 After 2016-08-02 2016-09-16 46
STTD 2017 Before 2017-07-14 2017-08-28 46
STTD 2017 During 2017-08-29 2017-09-01 4
STTD 2017 After 2017-09-20 2017-09-25 6
STTD 2018 Before 2018-07-20 2018-08-27 39
STTD 2018 During 2018-08-28 2018-09-26 30
STTD 2018 After 2018-09-27 2018-10-15 19
STTD 2019 Before 2019-07-26 2019-08-25 31
STTD 2019 During 2019-08-26 2019-09-21 27
STTD 2019 After 2019-09-22 2019-11-06 46
LIB 2015 Before 2015-07-06 2015-08-20 46
LIB 2015 During 2015-08-21 2015-10-01 38
LIB 2015 After 2015-10-02 2015-11-16 46
LIB 2016 Before 2016-05-29 2016-07-13 46
LIB 2016 During 2016-07-14 2016-08-01 19
LIB 2016 After 2016-08-02 2016-09-16 46
LIB 2017 Before 2017-07-14 2017-08-28 46
LIB 2017 During 2017-08-29 2017-09-18 21
LIB 2017 After 2017-09-19 2017-11-03 46
LIB 2018 Before 2018-08-14 2018-08-27 14
LIB 2018 During 2018-08-28 2018-09-26 30
LIB 2018 After 2018-09-27 2018-11-11 46
LIB 2019 Before 2019-07-11 2019-07-30 20
LIB 2019 During 2019-09-05 2019-09-21 17
LIB 2019 After 2019-09-22 2019-11-06 39
RVB 2015 Before 2015-07-06 2015-08-20 46
RVB 2015 During 2015-08-21 2015-10-01 42
RVB 2015 After 2015-10-02 2015-11-16 46
RVB 2016 Before 2016-05-29 2016-07-13 39
RVB 2016 During 2016-07-14 2016-08-01 19
RVB 2016 After 2016-08-02 2016-09-16 37
RVB 2017 Before 2017-07-14 2017-08-28 46
RVB 2017 During 2017-08-29 2017-09-18 21
RVB 2017 After 2017-09-19 2017-11-03 46
RVB 2018 Before 2018-07-13 2018-08-27 46
RVB 2018 During 2018-08-28 2018-09-26 30
RVB 2018 After 2018-09-27 2018-11-11 46
RVB 2019 Before 2019-07-11 2019-08-25 46
RVB 2019 During 2019-08-26 2019-09-21 27
RVB 2019 After 2019-09-22 2019-11-06 46

Before we run our models, we’ll need to remove the NA values in the Flow variable. When we remove all records where there are missing daily average flow values, there are even more Station-Year-FlowPeriod combinations with low sample counts within 2015-2019.

df_chla_wt_q <- df_chla_wt_c2 %>% drop_na(Flow)

df_chla_wt_q %>% 
  drop_na(Flow) %>% 
  summarize(
    min_date = min(Date),
    max_date = max(Date),
    num_samples = n(),
    .by = c(StationCode, Year, FlowActionPeriod)
  ) %>% 
  arrange(StationCode, Year, FlowActionPeriod) %>% 
  kable()
StationCode Year FlowActionPeriod min_date max_date num_samples
RCS 2015 Before 2015-08-17 2015-08-20 4
RCS 2015 During 2015-08-21 2015-09-28 39
RCS 2016 Before 2016-06-23 2016-07-13 8
RCS 2016 During 2016-07-14 2016-08-01 19
RCS 2016 After 2016-08-02 2016-09-16 46
RCS 2017 Before 2017-07-14 2017-08-28 46
RCS 2017 During 2017-08-29 2017-09-18 21
RCS 2017 After 2017-09-19 2017-11-03 36
RCS 2018 Before 2018-07-13 2018-08-27 46
RCS 2018 During 2018-08-28 2018-09-26 30
RCS 2018 After 2018-09-27 2018-11-11 46
RCS 2019 Before 2019-07-11 2019-08-25 46
RCS 2019 During 2019-08-26 2019-09-21 27
RCS 2019 After 2019-09-22 2019-11-06 46
RD22 2015 Before 2015-07-24 2015-08-20 28
RD22 2015 During 2015-08-21 2015-10-01 42
RD22 2015 After 2015-10-02 2015-11-10 40
RD22 2016 Before 2016-06-23 2016-07-13 21
RD22 2016 During 2016-07-14 2016-08-01 19
RD22 2016 After 2016-08-02 2016-09-16 46
RD22 2017 Before 2017-07-14 2017-08-28 46
RD22 2017 During 2017-08-29 2017-09-18 21
RD22 2017 After 2017-09-19 2017-11-03 46
RD22 2018 Before 2018-07-13 2018-08-27 46
RD22 2018 During 2018-08-28 2018-09-26 30
RD22 2018 After 2018-09-27 2018-11-11 46
RD22 2019 Before 2019-07-11 2019-08-25 46
RD22 2019 During 2019-08-26 2019-09-21 27
RD22 2019 After 2019-09-22 2019-11-06 46
I80 2015 Before 2015-07-24 2015-08-07 15
I80 2015 During 2015-08-27 2015-10-01 36
I80 2015 After 2015-10-02 2015-11-10 40
I80 2016 Before 2016-06-23 2016-07-13 21
I80 2016 During 2016-07-14 2016-08-01 19
I80 2016 After 2016-08-02 2016-09-16 46
I80 2017 Before 2017-07-14 2017-08-28 46
I80 2017 During 2017-08-29 2017-09-18 21
I80 2017 After 2017-09-19 2017-11-03 46
I80 2018 Before 2018-07-13 2018-08-27 46
I80 2018 During 2018-08-28 2018-09-26 30
I80 2018 After 2018-09-27 2018-11-11 46
I80 2019 Before 2019-07-11 2019-08-25 46
I80 2019 During 2019-08-26 2019-09-21 27
I80 2019 After 2019-09-22 2019-11-06 46
LIS 2015 Before 2015-07-06 2015-08-20 46
LIS 2015 During 2015-08-21 2015-10-01 37
LIS 2015 After 2015-10-02 2015-11-16 46
LIS 2016 Before 2016-05-29 2016-07-13 46
LIS 2016 During 2016-07-14 2016-08-01 19
LIS 2016 After 2016-08-02 2016-09-16 43
LIS 2017 Before 2017-07-14 2017-08-28 46
LIS 2017 During 2017-08-29 2017-09-09 12
LIS 2017 After 2017-09-19 2017-11-03 46
LIS 2018 Before 2018-07-13 2018-08-27 46
LIS 2018 During 2018-08-28 2018-09-26 30
LIS 2018 After 2018-09-27 2018-11-11 42
LIS 2019 Before 2019-07-19 2019-08-25 34
LIS 2019 During 2019-08-26 2019-09-21 27
LIS 2019 After 2019-09-22 2019-11-06 46
STTD 2015 Before 2015-07-27 2015-08-20 25
STTD 2015 During 2015-08-21 2015-10-01 37
STTD 2015 After 2015-10-02 2015-11-16 46
STTD 2016 Before 2016-06-23 2016-07-13 21
STTD 2016 During 2016-07-14 2016-08-01 19
STTD 2016 After 2016-08-02 2016-09-16 43
STTD 2017 Before 2017-07-14 2017-08-28 46
STTD 2017 During 2017-08-29 2017-09-01 4
STTD 2017 After 2017-09-20 2017-09-25 6
STTD 2018 Before 2018-07-20 2018-08-27 39
STTD 2018 During 2018-08-28 2018-09-26 30
STTD 2018 After 2018-09-27 2018-10-11 15
STTD 2019 Before 2019-07-26 2019-08-25 27
STTD 2019 During 2019-08-26 2019-09-21 27
STTD 2019 After 2019-09-22 2019-11-06 46
LIB 2015 Before 2015-07-06 2015-08-20 46
LIB 2015 During 2015-08-21 2015-10-01 33
LIB 2015 After 2015-10-02 2015-11-16 46
LIB 2016 Before 2016-05-29 2016-07-13 36
LIB 2016 During 2016-07-14 2016-08-01 19
LIB 2016 After 2016-08-02 2016-09-16 46
LIB 2017 Before 2017-07-14 2017-08-28 46
LIB 2017 During 2017-08-29 2017-09-18 21
LIB 2017 After 2017-09-19 2017-11-02 29
LIB 2018 Before 2018-08-14 2018-08-27 9
LIB 2018 During 2018-08-28 2018-09-08 12
LIB 2018 After 2018-09-27 2018-10-31 15
LIB 2019 Before 2019-07-11 2019-07-30 18
LIB 2019 During 2019-09-07 2019-09-21 10
LIB 2019 After 2019-09-22 2019-10-25 3
RVB 2015 Before 2015-07-06 2015-08-20 46
RVB 2015 During 2015-08-21 2015-10-01 42
RVB 2015 After 2015-10-02 2015-11-16 46
RVB 2016 Before 2016-05-29 2016-07-13 39
RVB 2016 During 2016-07-14 2016-08-01 4
RVB 2016 After 2016-08-02 2016-09-16 37
RVB 2017 Before 2017-07-14 2017-08-28 46
RVB 2017 During 2017-08-29 2017-09-18 21
RVB 2017 After 2017-09-19 2017-11-03 46
RVB 2018 Before 2018-07-13 2018-08-27 46
RVB 2018 During 2018-08-28 2018-09-26 30
RVB 2018 After 2018-09-27 2018-11-11 46
RVB 2019 Before 2019-07-11 2019-08-25 39
RVB 2019 During 2019-08-26 2019-09-21 27
RVB 2019 After 2019-09-22 2019-11-06 36

Models

GAM Model 1: Categorical predictors

Plots

df_chla_wt_q %>%
  ggplot(aes(x = FlowActionPeriod, y = Chla, fill = FlowActionPeriod)) +
  geom_boxplot() +
  facet_grid(rows = vars(Year), cols = vars(StationCode)) +
  scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
  annotation_logticks(sides = "l") +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    legend.position = "top"
  )

Initial Model

We’ll try running a GAM including all two-way interactions of the categorical predictors - Year, Flow Action Period, Station - and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.

m_gam_cat <- gam(
  Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5), 
  data = df_chla_wt_q,
  method = "REML"
)

Lets look at the model summary and diagnostics:

summary(m_gam_cat)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, 
##     k = 5)
## 
## Parametric coefficients:
##                                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             8.92533    0.08173 109.206  < 2e-16 ***
## Year_fct2016                           -0.12304    0.09641  -1.276 0.201987    
## Year_fct2017                            0.43975    0.08230   5.343 9.72e-08 ***
## Year_fct2018                            0.15352    0.08074   1.901 0.057351 .  
## Year_fct2019                            0.06346    0.08131   0.781 0.435138    
## FlowActionPeriodDuring                 -0.01236    0.06708  -0.184 0.853795    
## FlowActionPeriodAfter                   0.58644    0.07712   7.605 3.67e-14 ***
## StationCodeRD22                         0.70108    0.08498   8.250 2.23e-16 ***
## StationCodeI80                          0.84468    0.08759   9.644  < 2e-16 ***
## StationCodeLIS                         -0.43029    0.08360  -5.147 2.80e-07 ***
## StationCodeSTTD                        -0.68280    0.08732  -7.820 7.00e-15 ***
## StationCodeLIB                         -1.25628    0.08535 -14.720  < 2e-16 ***
## StationCodeRVB                         -1.28890    0.08343 -15.449  < 2e-16 ***
## Year_fct2016:FlowActionPeriodDuring    -0.11791    0.07122  -1.656 0.097897 .  
## Year_fct2017:FlowActionPeriodDuring    -0.07107    0.05402  -1.316 0.188334    
## Year_fct2018:FlowActionPeriodDuring    -0.03788    0.05025  -0.754 0.451023    
## Year_fct2019:FlowActionPeriodDuring    -0.48857    0.05151  -9.485  < 2e-16 ***
## Year_fct2016:FlowActionPeriodAfter     -0.16200    0.07631  -2.123 0.033843 *  
## Year_fct2017:FlowActionPeriodAfter     -0.26360    0.04869  -5.414 6.59e-08 ***
## Year_fct2018:FlowActionPeriodAfter     -0.15729    0.04867  -3.232 0.001241 ** 
## Year_fct2019:FlowActionPeriodAfter      0.06520    0.04933   1.322 0.186312    
## Year_fct2016:StationCodeRD22           -0.27091    0.09536  -2.841 0.004527 ** 
## Year_fct2017:StationCodeRD22           -0.60584    0.09026  -6.713 2.23e-11 ***
## Year_fct2018:StationCodeRD22           -0.42354    0.08782  -4.823 1.48e-06 ***
## Year_fct2019:StationCodeRD22           -0.21418    0.08845  -2.421 0.015512 *  
## Year_fct2016:StationCodeI80             0.35595    0.09680   3.677 0.000240 ***
## Year_fct2017:StationCodeI80            -0.23252    0.09193  -2.529 0.011478 *  
## Year_fct2018:StationCodeI80            -0.49323    0.08953  -5.509 3.88e-08 ***
## Year_fct2019:StationCodeI80             0.45168    0.09016   5.010 5.73e-07 ***
## Year_fct2016:StationCodeLIS             0.75033    0.09371   8.007 1.60e-15 ***
## Year_fct2017:StationCodeLIS             0.23373    0.09054   2.581 0.009882 ** 
## Year_fct2018:StationCodeLIS            -0.08130    0.08738  -0.930 0.352235    
## Year_fct2019:StationCodeLIS             0.63379    0.08855   7.157 1.00e-12 ***
## Year_fct2016:StationCodeSTTD            0.68029    0.09599   7.087 1.66e-12 ***
## Year_fct2017:StationCodeSTTD           -0.06767    0.10052  -0.673 0.500848    
## Year_fct2018:StationCodeSTTD           -0.61835    0.09139  -6.766 1.55e-11 ***
## Year_fct2019:StationCodeSTTD           -0.82420    0.08995  -9.163  < 2e-16 ***
## Year_fct2016:StationCodeLIB             1.17255    0.09423  12.443  < 2e-16 ***
## Year_fct2017:StationCodeLIB            -0.27805    0.09103  -3.054 0.002272 ** 
## Year_fct2018:StationCodeLIB            -1.92894    0.10180 -18.949  < 2e-16 ***
## Year_fct2019:StationCodeLIB            -0.61755    0.10553  -5.852 5.32e-09 ***
## Year_fct2016:StationCodeRVB             0.35867    0.09694   3.700 0.000219 ***
## Year_fct2017:StationCodeRVB            -0.55956    0.08929  -6.267 4.14e-10 ***
## Year_fct2018:StationCodeRVB            -0.38660    0.08685  -4.452 8.80e-06 ***
## Year_fct2019:StationCodeRVB            -0.65441    0.08842  -7.401 1.70e-13 ***
## FlowActionPeriodDuring:StationCodeRD22 -0.41717    0.06339  -6.581 5.40e-11 ***
## FlowActionPeriodAfter:StationCodeRD22  -0.62430    0.05636 -11.077  < 2e-16 ***
## FlowActionPeriodDuring:StationCodeI80  -0.51777    0.06421  -8.063 1.02e-15 ***
## FlowActionPeriodAfter:StationCodeI80   -0.50146    0.05684  -8.822  < 2e-16 ***
## FlowActionPeriodDuring:StationCodeLIS   0.42332    0.06385   6.630 3.90e-11 ***
## FlowActionPeriodAfter:StationCodeLIS   -0.52929    0.05589  -9.470  < 2e-16 ***
## FlowActionPeriodDuring:StationCodeSTTD  1.12064    0.06775  16.541  < 2e-16 ***
## FlowActionPeriodAfter:StationCodeSTTD  -0.38714    0.06322  -6.123 1.02e-09 ***
## FlowActionPeriodDuring:StationCodeLIB   0.26405    0.06912   3.820 0.000136 ***
## FlowActionPeriodAfter:StationCodeLIB   -0.75378    0.06336 -11.897  < 2e-16 ***
## FlowActionPeriodDuring:StationCodeRVB   0.08567    0.06416   1.335 0.181922    
## FlowActionPeriodAfter:StationCodeRVB   -0.69367    0.05617 -12.350  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.772  3.974 27.07  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.864   Deviance explained = 86.6%
## -REML = 1561.7  Scale est. = 0.13479   n = 3478
appraise(m_gam_cat)

k.check(m_gam_cat)
##        k'     edf   k-index p-value
## s(DOY)  4 3.77179 0.9457395       0
draw(m_gam_cat, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gam_cat, pages = 1, all.terms = TRUE)

acf(residuals(m_gam_cat))

Box.test(residuals(m_gam_cat), lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuals(m_gam_cat)
## X-squared = 8609.3, df = 20, p-value < 2.2e-16

Model with AR() correlation structure

Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.

# Define model formula as an object
f_gam_cat <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5)")

# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat <- gamm(
  f_gam_cat, 
  data = df_chla_wt_q,
  method = "REML"
)

m_gamm_cat_ar1 <- gamm(
  f_gam_cat, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_cat_ar2 <- gamm(
  f_gam_cat, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_cat_ar3 <- gamm(
  f_gam_cat, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
  method = "REML"
)

# Compare models
anova(
  m_gamm_cat$lme, 
  m_gamm_cat_ar1$lme, 
  m_gamm_cat_ar2$lme, 
  m_gamm_cat_ar3$lme
)
##                    Model df       AIC       BIC    logLik   Test  L.Ratio
## m_gamm_cat$lme         1 60  3243.372  3611.616 -1561.686                
## m_gamm_cat_ar1$lme     2 61 -2167.534 -1793.153  1144.767 1 vs 2 5412.906
## m_gamm_cat_ar2$lme     3 62 -2189.412 -1808.893  1156.706 2 vs 3   23.878
## m_gamm_cat_ar3$lme     4 63 -2188.415 -1801.759  1157.208 3 vs 4    1.004
##                    p-value
## m_gamm_cat$lme            
## m_gamm_cat_ar1$lme  <.0001
## m_gamm_cat_ar2$lme  <.0001
## m_gamm_cat_ar3$lme  0.3164

It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.

AR(2) Model

# Pull out the residuals and the GAM model
resid_cat_ar2 <- residuals(m_gamm_cat_ar2$lme, type = "normalized")
m_gamm_cat_ar2_gam <- m_gamm_cat_ar2$gam

summary(m_gamm_cat_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, 
##     k = 5)
## 
## Parametric coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             9.092494   0.347301  26.180  < 2e-16
## Year_fct2016                            0.074100   0.454595   0.163 0.870526
## Year_fct2017                            0.351443   0.428540   0.820 0.412220
## Year_fct2018                            0.220090   0.419794   0.524 0.600117
## Year_fct2019                           -0.004065   0.421429  -0.010 0.992305
## FlowActionPeriodDuring                  0.106229   0.094542   1.124 0.261253
## FlowActionPeriodAfter                   0.205035   0.143274   1.431 0.152500
## StationCodeRD22                         0.213335   0.427422   0.499 0.617727
## StationCodeI80                          0.541690   0.439999   1.231 0.218364
## StationCodeLIS                         -0.597751   0.419022  -1.427 0.153805
## StationCodeSTTD                        -0.309144   0.429632  -0.720 0.471848
## StationCodeLIB                         -1.392380   0.421407  -3.304 0.000963
## StationCodeRVB                         -1.450049   0.417001  -3.477 0.000513
## Year_fct2016:FlowActionPeriodDuring    -0.162149   0.090727  -1.787 0.073991
## Year_fct2017:FlowActionPeriodDuring     0.011836   0.089175   0.133 0.894417
## Year_fct2018:FlowActionPeriodDuring    -0.171794   0.089080  -1.929 0.053872
## Year_fct2019:FlowActionPeriodDuring     0.054852   0.089154   0.615 0.538428
## Year_fct2016:FlowActionPeriodAfter     -0.216376   0.127679  -1.695 0.090227
## Year_fct2017:FlowActionPeriodAfter     -0.082755   0.125185  -0.661 0.508619
## Year_fct2018:FlowActionPeriodAfter     -0.235315   0.125043  -1.882 0.059938
## Year_fct2019:FlowActionPeriodAfter      0.127736   0.124974   1.022 0.306806
## Year_fct2016:StationCodeRD22           -0.409097   0.576632  -0.709 0.478087
## Year_fct2017:StationCodeRD22           -0.510460   0.546951  -0.933 0.350740
## Year_fct2018:StationCodeRD22           -0.375917   0.536410  -0.701 0.483474
## Year_fct2019:StationCodeRD22           -0.125667   0.538556  -0.233 0.815511
## Year_fct2016:StationCodeI80             0.159128   0.585194   0.272 0.785697
## Year_fct2017:StationCodeI80            -0.229520   0.556103  -0.413 0.679830
## Year_fct2018:StationCodeI80            -0.577356   0.545768  -1.058 0.290186
## Year_fct2019:StationCodeI80             0.412105   0.547880   0.752 0.451994
## Year_fct2016:StationCodeLIS             0.766734   0.561958   1.364 0.172532
## Year_fct2017:StationCodeLIS             0.178449   0.544443   0.328 0.743110
## Year_fct2018:StationCodeLIS            -0.130811   0.531634  -0.246 0.805655
## Year_fct2019:StationCodeLIS             0.384296   0.537308   0.715 0.474519
## Year_fct2016:StationCodeSTTD            0.494794   0.579473   0.854 0.393237
## Year_fct2017:StationCodeSTTD           -0.666965   0.584959  -1.140 0.254287
## Year_fct2018:StationCodeSTTD           -0.891467   0.556482  -1.602 0.109255
## Year_fct2019:StationCodeSTTD           -1.068222   0.547691  -1.950 0.051209
## Year_fct2016:StationCodeLIB             1.032582   0.565695   1.825 0.068038
## Year_fct2017:StationCodeLIB            -0.298959   0.549421  -0.544 0.586385
## Year_fct2018:StationCodeLIB            -2.012752   0.597348  -3.369 0.000761
## Year_fct2019:StationCodeLIB            -0.761233   0.605732  -1.257 0.208942
## Year_fct2016:StationCodeRVB             0.385100   0.575289   0.669 0.503284
## Year_fct2017:StationCodeRVB            -0.578913   0.539279  -1.073 0.283125
## Year_fct2018:StationCodeRVB            -0.447074   0.528583  -0.846 0.397726
## Year_fct2019:StationCodeRVB            -0.854363   0.537540  -1.589 0.112064
## FlowActionPeriodDuring:StationCodeRD22 -0.067302   0.105406  -0.638 0.523192
## FlowActionPeriodAfter:StationCodeRD22   0.031609   0.148699   0.213 0.831675
## FlowActionPeriodDuring:StationCodeI80  -0.058480   0.105829  -0.553 0.580579
## FlowActionPeriodAfter:StationCodeI80   -0.002189   0.149048  -0.015 0.988284
## FlowActionPeriodDuring:StationCodeLIS  -0.087621   0.105289  -0.832 0.405357
## FlowActionPeriodAfter:StationCodeLIS    0.013527   0.148521   0.091 0.927435
## FlowActionPeriodDuring:StationCodeSTTD  0.056894   0.105964   0.537 0.591360
## FlowActionPeriodAfter:StationCodeSTTD  -0.493124   0.150394  -3.279 0.001053
## FlowActionPeriodDuring:StationCodeLIB  -0.098209   0.106665  -0.921 0.357264
## FlowActionPeriodAfter:StationCodeLIB   -0.224653   0.150868  -1.489 0.136561
## FlowActionPeriodDuring:StationCodeRVB  -0.067906   0.105332  -0.645 0.519175
## FlowActionPeriodAfter:StationCodeRVB   -0.174298   0.148565  -1.173 0.240792
##                                           
## (Intercept)                            ***
## Year_fct2016                              
## Year_fct2017                              
## Year_fct2018                              
## Year_fct2019                              
## FlowActionPeriodDuring                    
## FlowActionPeriodAfter                     
## StationCodeRD22                           
## StationCodeI80                            
## StationCodeLIS                            
## StationCodeSTTD                           
## StationCodeLIB                         ***
## StationCodeRVB                         ***
## Year_fct2016:FlowActionPeriodDuring    .  
## Year_fct2017:FlowActionPeriodDuring       
## Year_fct2018:FlowActionPeriodDuring    .  
## Year_fct2019:FlowActionPeriodDuring       
## Year_fct2016:FlowActionPeriodAfter     .  
## Year_fct2017:FlowActionPeriodAfter        
## Year_fct2018:FlowActionPeriodAfter     .  
## Year_fct2019:FlowActionPeriodAfter        
## Year_fct2016:StationCodeRD22              
## Year_fct2017:StationCodeRD22              
## Year_fct2018:StationCodeRD22              
## Year_fct2019:StationCodeRD22              
## Year_fct2016:StationCodeI80               
## Year_fct2017:StationCodeI80               
## Year_fct2018:StationCodeI80               
## Year_fct2019:StationCodeI80               
## Year_fct2016:StationCodeLIS               
## Year_fct2017:StationCodeLIS               
## Year_fct2018:StationCodeLIS               
## Year_fct2019:StationCodeLIS               
## Year_fct2016:StationCodeSTTD              
## Year_fct2017:StationCodeSTTD              
## Year_fct2018:StationCodeSTTD              
## Year_fct2019:StationCodeSTTD           .  
## Year_fct2016:StationCodeLIB            .  
## Year_fct2017:StationCodeLIB               
## Year_fct2018:StationCodeLIB            ***
## Year_fct2019:StationCodeLIB               
## Year_fct2016:StationCodeRVB               
## Year_fct2017:StationCodeRVB               
## Year_fct2018:StationCodeRVB               
## Year_fct2019:StationCodeRVB               
## FlowActionPeriodDuring:StationCodeRD22    
## FlowActionPeriodAfter:StationCodeRD22     
## FlowActionPeriodDuring:StationCodeI80     
## FlowActionPeriodAfter:StationCodeI80      
## FlowActionPeriodDuring:StationCodeLIS     
## FlowActionPeriodAfter:StationCodeLIS      
## FlowActionPeriodDuring:StationCodeSTTD    
## FlowActionPeriodAfter:StationCodeSTTD  ** 
## FlowActionPeriodDuring:StationCodeLIB     
## FlowActionPeriodAfter:StationCodeLIB      
## FlowActionPeriodDuring:StationCodeRVB     
## FlowActionPeriodAfter:StationCodeRVB      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.734  3.734 11.49  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.785   
##   Scale est. = 0.27538   n = 3478
appraise(m_gamm_cat_ar2_gam)

k.check(m_gamm_cat_ar2_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.734167 0.7698425       0
draw(m_gamm_cat_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_cat_ar2_gam, pages = 1, all.terms = TRUE)

acf(resid_cat_ar2)

Box.test(resid_cat_ar2, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_cat_ar2
## X-squared = 32.886, df = 20, p-value = 0.03473

The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.

# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, 
##     k = 5)
## 
## Parametric Terms:
##                              df     F  p-value
## Year_fct                      4 0.345   0.8477
## FlowActionPeriod              2 1.070   0.3432
## StationCode                   6 9.384 3.08e-10
## Year_fct:FlowActionPeriod     8 1.991   0.0438
## Year_fct:StationCode         24 2.732 1.16e-05
## FlowActionPeriod:StationCode 12 4.817 6.68e-08
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value
## s(DOY) 3.734  3.734 11.49  <2e-16

All the interaction terms are significant in this model, so we’ll use this one for the model selection process.

GAM Model 2: Categorical predictors adding water temperature

Initial Model

We’ll try running the categorical model adding water temperature to see if it improves the predictive power of the model.

m_gam_cat_wt <- gam(
  Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + s(DOY, k = 5), 
  data = df_chla_wt_q,
  method = "REML"
)

Lets look at the model summary and diagnostics:

summary(m_gam_cat_wt)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + 
##     s(DOY, k = 5)
## 
## Parametric coefficients:
##                                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             9.4732371  0.4077027  23.236  < 2e-16
## Year_fct2016                            0.1185304  0.3872590   0.306 0.759566
## Year_fct2017                            0.1743035  0.3365746   0.518 0.604579
## Year_fct2018                            0.3038120  0.2848691   1.066 0.286275
## Year_fct2019                            0.0558398  0.2655069   0.210 0.833435
## FlowActionPeriodDuring                 -0.8448474  0.3953310  -2.137 0.032664
## FlowActionPeriodAfter                   0.1409841  0.3375317   0.418 0.676199
## StationCodeRD22                        -0.0797142  0.3149937  -0.253 0.800233
## StationCodeI80                          1.0484912  0.3010805   3.482 0.000503
## StationCodeLIS                         -3.0959341  0.2998499 -10.325  < 2e-16
## StationCodeSTTD                        -2.1712782  0.3221696  -6.740 1.86e-11
## StationCodeLIB                         -2.4195633  0.3868708  -6.254 4.49e-10
## StationCodeRVB                         -2.5192484  0.3578309  -7.040 2.31e-12
## WaterTemp                              -0.0172525  0.0159383  -1.082 0.279127
## Year_fct2016:FlowActionPeriodDuring    -0.1478412  0.0717935  -2.059 0.039545
## Year_fct2017:FlowActionPeriodDuring    -0.1142243  0.0559969  -2.040 0.041444
## Year_fct2018:FlowActionPeriodDuring     0.0081453  0.0566778   0.144 0.885736
## Year_fct2019:FlowActionPeriodDuring    -0.4834680  0.0537729  -8.991  < 2e-16
## Year_fct2016:FlowActionPeriodAfter     -0.0463716  0.0939445  -0.494 0.621616
## Year_fct2017:FlowActionPeriodAfter     -0.1164966  0.0985549  -1.182 0.237269
## Year_fct2018:FlowActionPeriodAfter     -0.1779551  0.0818456  -2.174 0.029753
## Year_fct2019:FlowActionPeriodAfter      0.1434797  0.0868389   1.652 0.098575
## Year_fct2016:StationCodeRD22           -0.4030850  0.1061711  -3.797 0.000149
## Year_fct2017:StationCodeRD22           -0.6258445  0.0886093  -7.063 1.97e-12
## Year_fct2018:StationCodeRD22           -0.4032341  0.0861352  -4.681 2.96e-06
## Year_fct2019:StationCodeRD22           -0.2105029  0.0863540  -2.438 0.014833
## Year_fct2016:StationCodeI80             0.3764403  0.1078367   3.491 0.000488
## Year_fct2017:StationCodeI80            -0.2152108  0.0907540  -2.371 0.017778
## Year_fct2018:StationCodeI80            -0.5036823  0.0880693  -5.719 1.16e-08
## Year_fct2019:StationCodeI80             0.4707630  0.0885154   5.318 1.11e-07
## Year_fct2016:StationCodeLIS             0.3980072  0.1037768   3.835 0.000128
## Year_fct2017:StationCodeLIS             0.1331875  0.0890675   1.495 0.134915
## Year_fct2018:StationCodeLIS            -0.0665137  0.0856235  -0.777 0.437321
## Year_fct2019:StationCodeLIS             0.5866469  0.0867510   6.762 1.59e-11
## Year_fct2016:StationCodeSTTD            0.4647981  0.1075972   4.320 1.61e-05
## Year_fct2017:StationCodeSTTD           -0.1040733  0.0993597  -1.047 0.294971
## Year_fct2018:StationCodeSTTD           -0.5941584  0.0894294  -6.644 3.54e-11
## Year_fct2019:StationCodeSTTD           -0.8406972  0.0880556  -9.547  < 2e-16
## Year_fct2016:StationCodeLIB             1.0594550  0.1104642   9.591  < 2e-16
## Year_fct2017:StationCodeLIB            -0.2288879  0.0930776  -2.459 0.013978
## Year_fct2018:StationCodeLIB            -1.8694287  0.1011190 -18.487  < 2e-16
## Year_fct2019:StationCodeLIB            -0.5687776  0.1049190  -5.421 6.34e-08
## Year_fct2016:StationCodeRVB             0.2273758  0.1104979   2.058 0.039691
## Year_fct2017:StationCodeRVB            -0.4993343  0.0911536  -5.478 4.62e-08
## Year_fct2018:StationCodeRVB            -0.3330771  0.0868161  -3.837 0.000127
## Year_fct2019:StationCodeRVB            -0.5914447  0.0887079  -6.667 3.03e-11
## Year_fct2016:WaterTemp                 -0.0047186  0.0147895  -0.319 0.749709
## Year_fct2017:WaterTemp                  0.0102785  0.0131868   0.779 0.435767
## Year_fct2018:WaterTemp                 -0.0074263  0.0113164  -0.656 0.511709
## Year_fct2019:WaterTemp                 -0.0009458  0.0102805  -0.092 0.926706
## FlowActionPeriodDuring:StationCodeRD22 -0.3367835  0.0692219  -4.865 1.19e-06
## FlowActionPeriodAfter:StationCodeRD22  -0.3952011  0.1035785  -3.815 0.000138
## FlowActionPeriodDuring:StationCodeI80  -0.5058395  0.0700402  -7.222 6.28e-13
## FlowActionPeriodAfter:StationCodeI80   -0.5403872  0.1029498  -5.249 1.62e-07
## FlowActionPeriodDuring:StationCodeLIS   0.6066616  0.0696598   8.709  < 2e-16
## FlowActionPeriodAfter:StationCodeLIS    0.1670770  0.0997364   1.675 0.093989
## FlowActionPeriodDuring:StationCodeSTTD  1.2258165  0.0734088  16.699  < 2e-16
## FlowActionPeriodAfter:StationCodeSTTD   0.0006864  0.1043389   0.007 0.994751
## FlowActionPeriodDuring:StationCodeLIB   0.4149955  0.0840477   4.938 8.29e-07
## FlowActionPeriodAfter:StationCodeLIB   -0.4732155  0.1097713  -4.311 1.67e-05
## FlowActionPeriodDuring:StationCodeRVB   0.2264187  0.0775096   2.921 0.003510
## FlowActionPeriodAfter:StationCodeRVB   -0.3932928  0.1064178  -3.696 0.000223
## FlowActionPeriodDuring:WaterTemp        0.0304368  0.0159083   1.913 0.055797
## FlowActionPeriodAfter:WaterTemp         0.0026541  0.0125480   0.212 0.832498
## StationCodeRD22:WaterTemp               0.0307775  0.0121747   2.528 0.011517
## StationCodeI80:WaterTemp               -0.0098427  0.0115431  -0.853 0.393890
## StationCodeLIS:WaterTemp                0.1106623  0.0116009   9.539  < 2e-16
## StationCodeSTTD:WaterTemp               0.0615658  0.0127363   4.834 1.40e-06
## StationCodeLIB:WaterTemp                0.0474936  0.0159554   2.977 0.002935
## StationCodeRVB:WaterTemp                0.0499363  0.0143724   3.474 0.000518
##                                           
## (Intercept)                            ***
## Year_fct2016                              
## Year_fct2017                              
## Year_fct2018                              
## Year_fct2019                              
## FlowActionPeriodDuring                 *  
## FlowActionPeriodAfter                     
## StationCodeRD22                           
## StationCodeI80                         ***
## StationCodeLIS                         ***
## StationCodeSTTD                        ***
## StationCodeLIB                         ***
## StationCodeRVB                         ***
## WaterTemp                                 
## Year_fct2016:FlowActionPeriodDuring    *  
## Year_fct2017:FlowActionPeriodDuring    *  
## Year_fct2018:FlowActionPeriodDuring       
## Year_fct2019:FlowActionPeriodDuring    ***
## Year_fct2016:FlowActionPeriodAfter        
## Year_fct2017:FlowActionPeriodAfter        
## Year_fct2018:FlowActionPeriodAfter     *  
## Year_fct2019:FlowActionPeriodAfter     .  
## Year_fct2016:StationCodeRD22           ***
## Year_fct2017:StationCodeRD22           ***
## Year_fct2018:StationCodeRD22           ***
## Year_fct2019:StationCodeRD22           *  
## Year_fct2016:StationCodeI80            ***
## Year_fct2017:StationCodeI80            *  
## Year_fct2018:StationCodeI80            ***
## Year_fct2019:StationCodeI80            ***
## Year_fct2016:StationCodeLIS            ***
## Year_fct2017:StationCodeLIS               
## Year_fct2018:StationCodeLIS               
## Year_fct2019:StationCodeLIS            ***
## Year_fct2016:StationCodeSTTD           ***
## Year_fct2017:StationCodeSTTD              
## Year_fct2018:StationCodeSTTD           ***
## Year_fct2019:StationCodeSTTD           ***
## Year_fct2016:StationCodeLIB            ***
## Year_fct2017:StationCodeLIB            *  
## Year_fct2018:StationCodeLIB            ***
## Year_fct2019:StationCodeLIB            ***
## Year_fct2016:StationCodeRVB            *  
## Year_fct2017:StationCodeRVB            ***
## Year_fct2018:StationCodeRVB            ***
## Year_fct2019:StationCodeRVB            ***
## Year_fct2016:WaterTemp                    
## Year_fct2017:WaterTemp                    
## Year_fct2018:WaterTemp                    
## Year_fct2019:WaterTemp                    
## FlowActionPeriodDuring:StationCodeRD22 ***
## FlowActionPeriodAfter:StationCodeRD22  ***
## FlowActionPeriodDuring:StationCodeI80  ***
## FlowActionPeriodAfter:StationCodeI80   ***
## FlowActionPeriodDuring:StationCodeLIS  ***
## FlowActionPeriodAfter:StationCodeLIS   .  
## FlowActionPeriodDuring:StationCodeSTTD ***
## FlowActionPeriodAfter:StationCodeSTTD     
## FlowActionPeriodDuring:StationCodeLIB  ***
## FlowActionPeriodAfter:StationCodeLIB   ***
## FlowActionPeriodDuring:StationCodeRVB  ** 
## FlowActionPeriodAfter:StationCodeRVB   ***
## FlowActionPeriodDuring:WaterTemp       .  
## FlowActionPeriodAfter:WaterTemp           
## StationCodeRD22:WaterTemp              *  
## StationCodeI80:WaterTemp                  
## StationCodeLIS:WaterTemp               ***
## StationCodeSTTD:WaterTemp              ***
## StationCodeLIB:WaterTemp               ** 
## StationCodeRVB:WaterTemp               ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.731  3.966 23.94  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.87   Deviance explained = 87.3%
## -REML = 1517.2  Scale est. = 0.12822   n = 3478
appraise(m_gam_cat_wt)

k.check(m_gam_cat_wt)
##        k'      edf   k-index p-value
## s(DOY)  4 3.730833 0.9463198       0
draw(m_gam_cat_wt, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gam_cat_wt, pages = 1, all.terms = TRUE)

acf(residuals(m_gam_cat_wt))

Box.test(residuals(m_gam_cat_wt), lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuals(m_gam_cat_wt)
## X-squared = 8031.8, df = 20, p-value < 2.2e-16

Model with AR() correlation structure

Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.

# Define model formula as an object
f_gam_cat_wt <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + s(DOY, k = 5)")

# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat_wt <- gamm(
  f_gam_cat_wt, 
  data = df_chla_wt_q,
  method = "REML"
)

m_gamm_cat_wt_ar1 <- gamm(
  f_gam_cat_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_cat_wt_ar2 <- gamm(
  f_gam_cat_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_cat_wt_ar3 <- gamm(
  f_gam_cat_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
  method = "REML"
)

# Compare models
anova(
  m_gamm_cat_wt$lme, 
  m_gamm_cat_wt_ar1$lme, 
  m_gamm_cat_wt_ar2$lme, 
  m_gamm_cat_wt_ar3$lme
)
##                       Model df       AIC       BIC    logLik   Test  L.Ratio
## m_gamm_cat_wt$lme         1 73  3180.470  3628.222 -1517.235                
## m_gamm_cat_wt_ar1$lme     2 74 -2135.571 -1681.686  1141.786 1 vs 2 5318.041
## m_gamm_cat_wt_ar2$lme     3 75 -2160.486 -1700.467  1155.243 2 vs 3   26.915
## m_gamm_cat_wt_ar3$lme     4 76 -2159.702 -1693.550  1155.851 3 vs 4    1.216
##                       p-value
## m_gamm_cat_wt$lme            
## m_gamm_cat_wt_ar1$lme  <.0001
## m_gamm_cat_wt_ar2$lme  <.0001
## m_gamm_cat_wt_ar3$lme  0.2701

It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.

AR(2) Model

# Pull out the residuals and the GAM model
resid_cat_wt_ar2 <- residuals(m_gamm_cat_wt_ar2$lme, type = "normalized")
m_gamm_cat_wt_ar2_gam <- m_gamm_cat_wt_ar2$gam

summary(m_gamm_cat_wt_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + 
##     s(DOY, k = 5)
## 
## Parametric coefficients:
##                                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             9.1445090  0.6901525  13.250  < 2e-16
## Year_fct2016                            1.3780976  0.6921213   1.991  0.04655
## Year_fct2017                            0.0579152  0.6440942   0.090  0.92836
## Year_fct2018                           -0.3370015  0.6545952  -0.515  0.60671
## Year_fct2019                           -0.3900285  0.6187327  -0.630  0.52850
## FlowActionPeriodDuring                 -0.4663337  0.4240927  -1.100  0.27158
## FlowActionPeriodAfter                  -0.3836823  0.3632026  -1.056  0.29087
## StationCodeRD22                        -1.1868159  0.6842598  -1.734  0.08293
## StationCodeI80                          0.9164038  0.6678780   1.372  0.17012
## StationCodeLIS                         -1.6220417  0.6764787  -2.398  0.01655
## StationCodeSTTD                        -0.5586425  0.6987943  -0.799  0.42409
## StationCodeLIB                         -0.2983707  0.7674643  -0.389  0.69747
## StationCodeRVB                         -1.9098221  0.8359931  -2.284  0.02240
## WaterTemp                               0.0009966  0.0243656   0.041  0.96738
## Year_fct2016:FlowActionPeriodDuring    -0.1572496  0.0926170  -1.698  0.08963
## Year_fct2017:FlowActionPeriodDuring    -0.0258551  0.0919774  -0.281  0.77865
## Year_fct2018:FlowActionPeriodDuring    -0.1036021  0.0974311  -1.063  0.28770
## Year_fct2019:FlowActionPeriodDuring     0.0270928  0.0904924   0.299  0.76466
## Year_fct2016:FlowActionPeriodAfter     -0.1711207  0.1383100  -1.237  0.21609
## Year_fct2017:FlowActionPeriodAfter     -0.0309569  0.1287914  -0.240  0.81006
## Year_fct2018:FlowActionPeriodAfter     -0.1354965  0.1343868  -1.008  0.31340
## Year_fct2019:FlowActionPeriodAfter      0.1469567  0.1263950   1.163  0.24504
## Year_fct2016:StationCodeRD22           -0.6182502  0.6279173  -0.985  0.32489
## Year_fct2017:StationCodeRD22           -0.5837985  0.5953210  -0.981  0.32684
## Year_fct2018:StationCodeRD22           -0.3239681  0.5851972  -0.554  0.57989
## Year_fct2019:StationCodeRD22           -0.1547261  0.5869605  -0.264  0.79210
## Year_fct2016:StationCodeI80             0.2258097  0.6387648   0.354  0.72373
## Year_fct2017:StationCodeI80            -0.1772145  0.6064542  -0.292  0.77014
## Year_fct2018:StationCodeI80            -0.5434617  0.5963020  -0.911  0.36216
## Year_fct2019:StationCodeI80             0.4333058  0.5981682   0.724  0.46888
## Year_fct2016:StationCodeLIS             0.5871339  0.6122516   0.959  0.33764
## Year_fct2017:StationCodeLIS             0.1547944  0.5925727   0.261  0.79394
## Year_fct2018:StationCodeLIS            -0.1245806  0.5799582  -0.215  0.82993
## Year_fct2019:StationCodeLIS             0.3524021  0.5857178   0.602  0.54744
## Year_fct2016:StationCodeSTTD            0.4354888  0.6320048   0.689  0.49083
## Year_fct2017:StationCodeSTTD           -0.6528634  0.6368032  -1.025  0.30533
## Year_fct2018:StationCodeSTTD           -0.8898961  0.6067842  -1.467  0.14258
## Year_fct2019:StationCodeSTTD           -1.0781623  0.5972945  -1.805  0.07115
## Year_fct2016:StationCodeLIB             0.9323841  0.6169532   1.511  0.13081
## Year_fct2017:StationCodeLIB            -0.2392731  0.5986330  -0.400  0.68940
## Year_fct2018:StationCodeLIB            -2.0226118  0.6499257  -3.112  0.00187
## Year_fct2019:StationCodeLIB            -0.7318547  0.6572217  -1.114  0.26555
## Year_fct2016:StationCodeRVB             0.2238684  0.6270456   0.357  0.72110
## Year_fct2017:StationCodeRVB            -0.5138895  0.5877741  -0.874  0.38202
## Year_fct2018:StationCodeRVB            -0.3928585  0.5778469  -0.680  0.49664
## Year_fct2019:StationCodeRVB            -0.7921245  0.5874493  -1.348  0.17762
## Year_fct2016:WaterTemp                 -0.0524611  0.0207223  -2.532  0.01140
## Year_fct2017:WaterTemp                  0.0126743  0.0191452   0.662  0.50801
## Year_fct2018:WaterTemp                  0.0245102  0.0205942   1.190  0.23407
## Year_fct2019:WaterTemp                  0.0193582  0.0181899   1.064  0.28730
## FlowActionPeriodDuring:StationCodeRD22 -0.0578678  0.1049823  -0.551  0.58152
## FlowActionPeriodAfter:StationCodeRD22   0.0980816  0.1507379   0.651  0.51530
## FlowActionPeriodDuring:StationCodeI80  -0.0674505  0.1057957  -0.638  0.52381
## FlowActionPeriodAfter:StationCodeI80   -0.0037150  0.1519709  -0.024  0.98050
## FlowActionPeriodDuring:StationCodeLIS  -0.0502238  0.1055792  -0.476  0.63432
## FlowActionPeriodAfter:StationCodeLIS    0.1157942  0.1515187   0.764  0.44479
## FlowActionPeriodDuring:StationCodeSTTD  0.0666638  0.1064582   0.626  0.53123
## FlowActionPeriodAfter:StationCodeSTTD  -0.4517271  0.1532547  -2.948  0.00322
## FlowActionPeriodDuring:StationCodeLIB  -0.0730926  0.1151386  -0.635  0.52559
## FlowActionPeriodAfter:StationCodeLIB   -0.2143517  0.1567062  -1.368  0.17145
## FlowActionPeriodDuring:StationCodeRVB  -0.0110737  0.1127300  -0.098  0.92175
## FlowActionPeriodAfter:StationCodeRVB   -0.0964396  0.1546565  -0.624  0.53295
## FlowActionPeriodDuring:WaterTemp        0.0231542  0.0168338   1.375  0.16908
## FlowActionPeriodAfter:WaterTemp         0.0210524  0.0143386   1.468  0.14213
## StationCodeRD22:WaterTemp               0.0653478  0.0222444   2.938  0.00333
## StationCodeI80:WaterTemp               -0.0191793  0.0207291  -0.925  0.35491
## StationCodeLIS:WaterTemp                0.0464759  0.0218848   2.124  0.03377
## StationCodeSTTD:WaterTemp               0.0103409  0.0232382   0.445  0.65635
## StationCodeLIB:WaterTemp               -0.0551245  0.0276723  -1.992  0.04645
## StationCodeRVB:WaterTemp                0.0192565  0.0319521   0.603  0.54677
##                                           
## (Intercept)                            ***
## Year_fct2016                           *  
## Year_fct2017                              
## Year_fct2018                              
## Year_fct2019                              
## FlowActionPeriodDuring                    
## FlowActionPeriodAfter                     
## StationCodeRD22                        .  
## StationCodeI80                            
## StationCodeLIS                         *  
## StationCodeSTTD                           
## StationCodeLIB                            
## StationCodeRVB                         *  
## WaterTemp                                 
## Year_fct2016:FlowActionPeriodDuring    .  
## Year_fct2017:FlowActionPeriodDuring       
## Year_fct2018:FlowActionPeriodDuring       
## Year_fct2019:FlowActionPeriodDuring       
## Year_fct2016:FlowActionPeriodAfter        
## Year_fct2017:FlowActionPeriodAfter        
## Year_fct2018:FlowActionPeriodAfter        
## Year_fct2019:FlowActionPeriodAfter        
## Year_fct2016:StationCodeRD22              
## Year_fct2017:StationCodeRD22              
## Year_fct2018:StationCodeRD22              
## Year_fct2019:StationCodeRD22              
## Year_fct2016:StationCodeI80               
## Year_fct2017:StationCodeI80               
## Year_fct2018:StationCodeI80               
## Year_fct2019:StationCodeI80               
## Year_fct2016:StationCodeLIS               
## Year_fct2017:StationCodeLIS               
## Year_fct2018:StationCodeLIS               
## Year_fct2019:StationCodeLIS               
## Year_fct2016:StationCodeSTTD              
## Year_fct2017:StationCodeSTTD              
## Year_fct2018:StationCodeSTTD              
## Year_fct2019:StationCodeSTTD           .  
## Year_fct2016:StationCodeLIB               
## Year_fct2017:StationCodeLIB               
## Year_fct2018:StationCodeLIB            ** 
## Year_fct2019:StationCodeLIB               
## Year_fct2016:StationCodeRVB               
## Year_fct2017:StationCodeRVB               
## Year_fct2018:StationCodeRVB               
## Year_fct2019:StationCodeRVB               
## Year_fct2016:WaterTemp                 *  
## Year_fct2017:WaterTemp                    
## Year_fct2018:WaterTemp                    
## Year_fct2019:WaterTemp                    
## FlowActionPeriodDuring:StationCodeRD22    
## FlowActionPeriodAfter:StationCodeRD22     
## FlowActionPeriodDuring:StationCodeI80     
## FlowActionPeriodAfter:StationCodeI80      
## FlowActionPeriodDuring:StationCodeLIS     
## FlowActionPeriodAfter:StationCodeLIS      
## FlowActionPeriodDuring:StationCodeSTTD    
## FlowActionPeriodAfter:StationCodeSTTD  ** 
## FlowActionPeriodDuring:StationCodeLIB     
## FlowActionPeriodAfter:StationCodeLIB      
## FlowActionPeriodDuring:StationCodeRVB     
## FlowActionPeriodAfter:StationCodeRVB      
## FlowActionPeriodDuring:WaterTemp          
## FlowActionPeriodAfter:WaterTemp           
## StationCodeRD22:WaterTemp              ** 
## StationCodeI80:WaterTemp                  
## StationCodeLIS:WaterTemp               *  
## StationCodeSTTD:WaterTemp                 
## StationCodeLIB:WaterTemp               *  
## StationCodeRVB:WaterTemp                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value    
## s(DOY) 3.669  3.669 8.849 5.92e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.774   
##   Scale est. = 0.30258   n = 3478
appraise(m_gamm_cat_wt_ar2_gam)

k.check(m_gamm_cat_wt_ar2_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.669285 0.7410218       0
draw(m_gamm_cat_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_cat_wt_ar2_gam, pages = 1, all.terms = TRUE)

acf(resid_cat_wt_ar2)

Box.test(resid_cat_wt_ar2, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_cat_wt_ar2
## X-squared = 30.859, df = 20, p-value = 0.05708

The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.

# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat_wt_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + 
##     s(DOY, k = 5)
## 
## Parametric Terms:
##                              df     F  p-value
## Year_fct                      4 2.412 0.047021
## FlowActionPeriod              2 0.740 0.477314
## StationCode                   6 5.258 2.11e-05
## WaterTemp                     1 0.002 0.967377
## Year_fct:FlowActionPeriod     8 1.159 0.320505
## Year_fct:StationCode         24 2.285 0.000357
## Year_fct:WaterTemp            4 4.446 0.001386
## FlowActionPeriod:StationCode 12 4.941 3.63e-08
## FlowActionPeriod:WaterTemp    2 1.414 0.243334
## StationCode:WaterTemp         6 6.427 9.53e-07
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value
## s(DOY) 3.669  3.669 8.849 5.92e-06

The Year:StationCode, Year:WaterTemp, FlowActionPeriod:StationCode, and StationCode:WaterTemp interactions were significant among the interaction terms of the model. Let’s remove the Year:FlowActionPeriod and FlowActionPeriod:WaterTemp interaction terms and check the model summary and diagnostics.

m_gamm_cat_wt_ar2b <- gamm(
  Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * FlowActionPeriod + s(DOY, k = 5), 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

# Pull out the residuals and the GAM model
resid_cat_wt_ar2b <- residuals(m_gamm_cat_wt_ar2b$lme, type = "normalized")
m_gamm_cat_wt_ar2b_gam <- m_gamm_cat_wt_ar2b$gam

summary(m_gamm_cat_wt_ar2b_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * 
##     FlowActionPeriod + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             8.782805   0.664490  13.217  < 2e-16
## Year_fct2016                            1.341143   0.707272   1.896  0.05802
## Year_fct2017                           -0.006022   0.652080  -0.009  0.99263
## Year_fct2018                           -0.528205   0.657513  -0.803  0.42184
## Year_fct2019                           -0.272579   0.628678  -0.434  0.66462
## StationCodeRD22                        -1.163183   0.707102  -1.645  0.10006
## StationCodeI80                          0.951126   0.691075   1.376  0.16882
## StationCodeLIS                         -1.559134   0.698191  -2.233  0.02561
## StationCodeSTTD                        -0.493857   0.720065  -0.686  0.49285
## StationCodeLIB                         -0.085665   0.784497  -0.109  0.91305
## StationCodeRVB                         -1.782352   0.855977  -2.082  0.03739
## WaterTemp                               0.016614   0.022904   0.725  0.46828
## FlowActionPeriodDuring                  0.052204   0.074987   0.696  0.48637
## FlowActionPeriodAfter                   0.099477   0.110745   0.898  0.36912
## Year_fct2016:StationCodeRD22           -0.613100   0.667140  -0.919  0.35816
## Year_fct2017:StationCodeRD22           -0.578281   0.633452  -0.913  0.36136
## Year_fct2018:StationCodeRD22           -0.315304   0.622783  -0.506  0.61269
## Year_fct2019:StationCodeRD22           -0.149630   0.624636  -0.240  0.81070
## Year_fct2016:StationCodeI80             0.223843   0.678341   0.330  0.74143
## Year_fct2017:StationCodeI80            -0.181161   0.644952  -0.281  0.77881
## Year_fct2018:StationCodeI80            -0.547381   0.634232  -0.863  0.38817
## Year_fct2019:StationCodeI80             0.423004   0.636191   0.665  0.50616
## Year_fct2016:StationCodeLIS             0.619479   0.650406   0.952  0.34094
## Year_fct2017:StationCodeLIS             0.170080   0.630490   0.270  0.78736
## Year_fct2018:StationCodeLIS            -0.116806   0.617195  -0.189  0.84991
## Year_fct2019:StationCodeLIS             0.349457   0.623223   0.561  0.57502
## Year_fct2016:StationCodeSTTD            0.445998   0.671203   0.664  0.50643
## Year_fct2017:StationCodeSTTD           -0.666675   0.676269  -0.986  0.32429
## Year_fct2018:StationCodeSTTD           -0.889459   0.645302  -1.378  0.16818
## Year_fct2019:StationCodeSTTD           -1.078389   0.635280  -1.698  0.08969
## Year_fct2016:StationCodeLIB             0.940927   0.655358   1.436  0.15117
## Year_fct2017:StationCodeLIB            -0.223119   0.636768  -0.350  0.72607
## Year_fct2018:StationCodeLIB            -2.019727   0.689265  -2.930  0.00341
## Year_fct2019:StationCodeLIB            -0.746492   0.696748  -1.071  0.28407
## Year_fct2016:StationCodeRVB             0.239251   0.665883   0.359  0.71939
## Year_fct2017:StationCodeRVB            -0.506117   0.625396  -0.809  0.41841
## Year_fct2018:StationCodeRVB            -0.384654   0.614812  -0.626  0.53159
## Year_fct2019:StationCodeRVB            -0.807469   0.624907  -1.292  0.19640
## Year_fct2016:WaterTemp                 -0.056611   0.020507  -2.761  0.00580
## Year_fct2017:WaterTemp                  0.013897   0.018782   0.740  0.45940
## Year_fct2018:WaterTemp                  0.029403   0.020156   1.459  0.14472
## Year_fct2019:WaterTemp                  0.016339   0.017903   0.913  0.36148
## StationCodeRD22:WaterTemp               0.064117   0.022486   2.851  0.00438
## StationCodeI80:WaterTemp               -0.020588   0.020915  -0.984  0.32500
## StationCodeLIS:WaterTemp                0.043491   0.022089   1.969  0.04904
## StationCodeSTTD:WaterTemp               0.008117   0.023400   0.347  0.72870
## StationCodeLIB:WaterTemp               -0.063750   0.027815  -2.292  0.02197
## StationCodeRVB:WaterTemp                0.015016   0.032392   0.464  0.64298
## StationCodeRD22:FlowActionPeriodDuring -0.054000   0.105338  -0.513  0.60824
## StationCodeI80:FlowActionPeriodDuring  -0.062165   0.106067  -0.586  0.55785
## StationCodeLIS:FlowActionPeriodDuring  -0.061674   0.105451  -0.585  0.55868
## StationCodeSTTD:FlowActionPeriodDuring  0.042061   0.105724   0.398  0.69078
## StationCodeLIB:FlowActionPeriodDuring  -0.133669   0.107077  -1.248  0.21199
## StationCodeRVB:FlowActionPeriodDuring  -0.063906   0.105456  -0.606  0.54456
## StationCodeRD22:FlowActionPeriodAfter   0.101187   0.150452   0.673  0.50128
## StationCodeI80:FlowActionPeriodAfter    0.005898   0.150783   0.039  0.96880
## StationCodeLIS:FlowActionPeriodAfter    0.107314   0.151165   0.710  0.47781
## StationCodeSTTD:FlowActionPeriodAfter  -0.476380   0.152352  -3.127  0.00178
## StationCodeLIB:FlowActionPeriodAfter   -0.265101   0.153142  -1.731  0.08353
## StationCodeRVB:FlowActionPeriodAfter   -0.142648   0.151288  -0.943  0.34580
##                                           
## (Intercept)                            ***
## Year_fct2016                           .  
## Year_fct2017                              
## Year_fct2018                              
## Year_fct2019                              
## StationCodeRD22                           
## StationCodeI80                            
## StationCodeLIS                         *  
## StationCodeSTTD                           
## StationCodeLIB                            
## StationCodeRVB                         *  
## WaterTemp                                 
## FlowActionPeriodDuring                    
## FlowActionPeriodAfter                     
## Year_fct2016:StationCodeRD22              
## Year_fct2017:StationCodeRD22              
## Year_fct2018:StationCodeRD22              
## Year_fct2019:StationCodeRD22              
## Year_fct2016:StationCodeI80               
## Year_fct2017:StationCodeI80               
## Year_fct2018:StationCodeI80               
## Year_fct2019:StationCodeI80               
## Year_fct2016:StationCodeLIS               
## Year_fct2017:StationCodeLIS               
## Year_fct2018:StationCodeLIS               
## Year_fct2019:StationCodeLIS               
## Year_fct2016:StationCodeSTTD              
## Year_fct2017:StationCodeSTTD              
## Year_fct2018:StationCodeSTTD              
## Year_fct2019:StationCodeSTTD           .  
## Year_fct2016:StationCodeLIB               
## Year_fct2017:StationCodeLIB               
## Year_fct2018:StationCodeLIB            ** 
## Year_fct2019:StationCodeLIB               
## Year_fct2016:StationCodeRVB               
## Year_fct2017:StationCodeRVB               
## Year_fct2018:StationCodeRVB               
## Year_fct2019:StationCodeRVB               
## Year_fct2016:WaterTemp                 ** 
## Year_fct2017:WaterTemp                    
## Year_fct2018:WaterTemp                    
## Year_fct2019:WaterTemp                    
## StationCodeRD22:WaterTemp              ** 
## StationCodeI80:WaterTemp                  
## StationCodeLIS:WaterTemp               *  
## StationCodeSTTD:WaterTemp                 
## StationCodeLIB:WaterTemp               *  
## StationCodeRVB:WaterTemp                  
## StationCodeRD22:FlowActionPeriodDuring    
## StationCodeI80:FlowActionPeriodDuring     
## StationCodeLIS:FlowActionPeriodDuring     
## StationCodeSTTD:FlowActionPeriodDuring    
## StationCodeLIB:FlowActionPeriodDuring     
## StationCodeRVB:FlowActionPeriodDuring     
## StationCodeRD22:FlowActionPeriodAfter     
## StationCodeI80:FlowActionPeriodAfter      
## StationCodeLIS:FlowActionPeriodAfter      
## StationCodeSTTD:FlowActionPeriodAfter  ** 
## StationCodeLIB:FlowActionPeriodAfter   .  
## StationCodeRVB:FlowActionPeriodAfter      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value    
## s(DOY) 3.656  3.656 7.364 1.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.762   
##   Scale est. = 0.32529   n = 3478
appraise(m_gamm_cat_wt_ar2b_gam)

k.check(m_gamm_cat_wt_ar2b_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.655508 0.7254327       0
draw(m_gamm_cat_wt_ar2b_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_cat_wt_ar2b_gam, pages = 1, all.terms = TRUE)

acf(resid_cat_wt_ar2b)

Box.test(resid_cat_wt_ar2b, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_cat_wt_ar2b
## X-squared = 31.046, df = 20, p-value = 0.05459
anova(m_gamm_cat_wt_ar2b_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * 
##     FlowActionPeriod + s(DOY, k = 5)
## 
## Parametric Terms:
##                              df     F  p-value
## Year_fct                      4 2.316 0.055024
## StationCode                   6 4.867 5.84e-05
## WaterTemp                     1 0.526 0.468278
## FlowActionPeriod              2 0.415 0.660676
## Year_fct:StationCode         24 2.036 0.002086
## Year_fct:WaterTemp            4 5.222 0.000342
## StationCode:WaterTemp         6 6.644 5.32e-07
## StationCode:FlowActionPeriod 12 5.012 2.55e-08
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value
## s(DOY) 3.656  3.656 7.364 1.34e-05

This model still looks pretty good. We’ll use this one for the model selection process.

GAM Model 3: Flow as continuous predictor

Plots

df_chla_wt_q %>% 
  ggplot(aes(x = Flow, y = Chla)) +
  geom_point() +
  geom_smooth(formula = "y ~ x") +
  facet_wrap(vars(StationCode), scales = "free") +
  theme_bw()

Let’s break these scatterplots apart by year to see how these patterns vary annually.

df_chla_wt_q %>% 
  ggplot(aes(x = Flow, y = Chla)) +
  geom_point() +
  geom_smooth(formula = "y ~ x") +
  facet_wrap(
    vars(StationCode, Year), 
    ncol = 5, 
    scales = "free", 
    labeller = labeller(.multi_line = FALSE)
  ) +
  theme_bw()

Initial Model

We’ll try running a GAM including all two-way interactions between Year, Daily Average Flow, and Station, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.

m_gam_q <- gam(
  Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5), 
  data = df_chla_wt_q,
  method = "REML"
)

Lets look at the model summary and diagnostics:

summary(m_gam_q)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   9.322e+00  6.781e-02 137.468  < 2e-16 ***
## Year_fct2016                 -6.414e-02  7.625e-02  -0.841 0.400279    
## Year_fct2017                  1.818e-01  7.628e-02   2.383 0.017221 *  
## Year_fct2018                  1.125e-01  7.096e-02   1.585 0.113065    
## Year_fct2019                  1.372e-02  7.116e-02   0.193 0.847183    
## Flow                         -1.275e-03  1.059e-04 -12.038  < 2e-16 ***
## StationCodeRD22               3.344e-01  7.804e-02   4.285 1.88e-05 ***
## StationCodeI80                5.559e-01  7.984e-02   6.963 3.98e-12 ***
## StationCodeLIS               -6.816e-01  7.585e-02  -8.987  < 2e-16 ***
## StationCodeSTTD              -6.447e-01  7.727e-02  -8.344  < 2e-16 ***
## StationCodeLIB               -1.640e+00  1.260e-01 -13.019  < 2e-16 ***
## StationCodeRVB               -1.908e+00  1.136e-01 -16.801  < 2e-16 ***
## Year_fct2016:Flow            -8.355e-05  3.159e-05  -2.645 0.008212 ** 
## Year_fct2017:Flow            -7.085e-05  2.628e-05  -2.696 0.007046 ** 
## Year_fct2018:Flow            -2.494e-05  2.636e-05  -0.946 0.344219    
## Year_fct2019:Flow            -5.410e-05  2.800e-05  -1.932 0.053386 .  
## Year_fct2016:StationCodeRD22 -4.235e-01  9.371e-02  -4.519 6.41e-06 ***
## Year_fct2017:StationCodeRD22 -5.668e-01  9.263e-02  -6.119 1.05e-09 ***
## Year_fct2018:StationCodeRD22 -4.231e-01  8.681e-02  -4.874 1.14e-06 ***
## Year_fct2019:StationCodeRD22 -1.754e-01  8.709e-02  -2.014 0.044067 *  
## Year_fct2016:StationCodeI80   2.386e-01  9.525e-02   2.505 0.012300 *  
## Year_fct2017:StationCodeI80  -2.309e-01  9.413e-02  -2.453 0.014198 *  
## Year_fct2018:StationCodeI80  -4.800e-01  8.848e-02  -5.425 6.21e-08 ***
## Year_fct2019:StationCodeI80   5.216e-01  8.876e-02   5.876 4.60e-09 ***
## Year_fct2016:StationCodeLIS   5.173e-01  9.177e-02   5.637 1.87e-08 ***
## Year_fct2017:StationCodeLIS   3.056e-01  9.172e-02   3.332 0.000870 ***
## Year_fct2018:StationCodeLIS  -1.288e-01  8.627e-02  -1.493 0.135495    
## Year_fct2019:StationCodeLIS   5.536e-01  8.765e-02   6.316 3.02e-10 ***
## Year_fct2016:StationCodeSTTD  3.119e-01  9.445e-02   3.302 0.000971 ***
## Year_fct2017:StationCodeSTTD -1.372e-01  9.943e-02  -1.380 0.167655    
## Year_fct2018:StationCodeSTTD -7.586e-01  9.080e-02  -8.354  < 2e-16 ***
## Year_fct2019:StationCodeSTTD -1.145e+00  8.934e-02 -12.815  < 2e-16 ***
## Year_fct2016:StationCodeLIB   8.329e-01  1.058e-01   7.871 4.67e-15 ***
## Year_fct2017:StationCodeLIB  -1.933e-01  1.063e-01  -1.818 0.069088 .  
## Year_fct2018:StationCodeLIB  -1.979e+00  1.223e-01 -16.186  < 2e-16 ***
## Year_fct2019:StationCodeLIB  -7.144e-01  1.241e-01  -5.757 9.31e-09 ***
## Year_fct2016:StationCodeRVB   6.335e-01  2.347e-01   2.699 0.006986 ** 
## Year_fct2017:StationCodeRVB  -2.198e-02  1.987e-01  -0.111 0.911935    
## Year_fct2018:StationCodeRVB  -4.140e-01  1.626e-01  -2.546 0.010955 *  
## Year_fct2019:StationCodeRVB  -4.845e-01  1.910e-01  -2.537 0.011238 *  
## Flow:StationCodeRD22         -1.057e-04  1.375e-04  -0.768 0.442446    
## Flow:StationCodeI80          -6.752e-04  1.387e-04  -4.867 1.18e-06 ***
## Flow:StationCodeLIS           1.541e-03  1.391e-04  11.079  < 2e-16 ***
## Flow:StationCodeSTTD          3.211e-03  1.406e-04  22.832  < 2e-16 ***
## Flow:StationCodeLIB           1.275e-03  1.220e-04  10.452  < 2e-16 ***
## Flow:StationCodeRVB           1.327e-03  1.039e-04  12.773  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.757  3.966 81.39  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.851   Deviance explained = 85.3%
## -REML = 1770.4  Scale est. = 0.14725   n = 3478
appraise(m_gam_q)

k.check(m_gam_q)
##        k'      edf   k-index p-value
## s(DOY)  4 3.756521 0.9228792       0
draw(m_gam_q, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gam_q, pages = 1, all.terms = TRUE)

acf(residuals(m_gam_q))

Box.test(residuals(m_gam_q), lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuals(m_gam_q)
## X-squared = 10096, df = 20, p-value < 2.2e-16

Model with AR() correlation structure

Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.

# Define model formula as an object
f_gam_q <- as.formula("Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)")

# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q <- gamm(
  f_gam_q, 
  data = df_chla_wt_q,
  method = "REML"
)

m_gamm_q_ar1 <- gamm(
  f_gam_q, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_ar2 <- gamm(
  f_gam_q, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_ar3 <- gamm(
  f_gam_q, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
  method = "REML"
)

# Compare models
anova(
  m_gamm_q$lme, 
  m_gamm_q_ar1$lme, 
  m_gamm_q_ar2$lme, 
  m_gamm_q_ar3$lme
)
##                  Model df       AIC       BIC    logLik   Test  L.Ratio p-value
## m_gamm_q$lme         1 49  3638.827  3939.716 -1770.413                        
## m_gamm_q_ar1$lme     2 50 -2092.217 -1785.187  1096.109 1 vs 2 5733.044  <.0001
## m_gamm_q_ar2$lme     3 51 -2111.366 -1798.195  1106.683 2 vs 3   21.149  <.0001
## m_gamm_q_ar3$lme     4 52 -2110.010 -1790.698  1107.005 3 vs 4    0.644  0.4224

It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.

AR(2) Model

# Pull out the residuals and the GAM model
resid_q_ar2 <- residuals(m_gamm_q_ar2$lme, type = "normalized")
m_gamm_q_ar2_gam <- m_gamm_q_ar2$gam

summary(m_gamm_q_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   9.178e+00  2.921e-01  31.417  < 2e-16 ***
## Year_fct2016                 -3.270e-02  3.761e-01  -0.087 0.930719    
## Year_fct2017                  3.329e-01  3.587e-01   0.928 0.353319    
## Year_fct2018                  1.306e-01  3.495e-01   0.374 0.708602    
## Year_fct2019                  8.786e-02  3.505e-01   0.251 0.802084    
## Flow                         -1.821e-04  1.772e-04  -1.028 0.304126    
## StationCodeRD22               3.170e-01  3.570e-01   0.888 0.374697    
## StationCodeI80                7.103e-01  3.671e-01   1.935 0.053071 .  
## StationCodeLIS               -5.832e-01  3.496e-01  -1.668 0.095356 .  
## StationCodeSTTD              -4.766e-01  3.576e-01  -1.333 0.182705    
## StationCodeLIB               -1.228e+00  3.569e-01  -3.440 0.000588 ***
## StationCodeRVB               -1.517e+00  3.519e-01  -4.312 1.66e-05 ***
## Year_fct2016:Flow            -5.265e-06  1.924e-05  -0.274 0.784344    
## Year_fct2017:Flow            -4.817e-06  1.964e-05  -0.245 0.806252    
## Year_fct2018:Flow            -3.172e-06  2.102e-05  -0.151 0.880067    
## Year_fct2019:Flow            -2.436e-05  1.961e-05  -1.242 0.214151    
## Year_fct2016:StationCodeRD22 -4.412e-01  4.823e-01  -0.915 0.360450    
## Year_fct2017:StationCodeRD22 -5.891e-01  4.591e-01  -1.283 0.199529    
## Year_fct2018:StationCodeRD22 -4.262e-01  4.486e-01  -0.950 0.342068    
## Year_fct2019:StationCodeRD22 -1.675e-01  4.503e-01  -0.372 0.709983    
## Year_fct2016:StationCodeI80   1.284e-01  4.896e-01   0.262 0.793087    
## Year_fct2017:StationCodeI80  -3.791e-01  4.669e-01  -0.812 0.416906    
## Year_fct2018:StationCodeI80  -6.219e-01  4.564e-01  -1.363 0.173092    
## Year_fct2019:StationCodeI80   3.949e-01  4.582e-01   0.862 0.388749    
## Year_fct2016:StationCodeLIS   6.885e-01  4.700e-01   1.465 0.143014    
## Year_fct2017:StationCodeLIS   1.647e-01  4.566e-01   0.361 0.718331    
## Year_fct2018:StationCodeLIS  -1.577e-01  4.445e-01  -0.355 0.722725    
## Year_fct2019:StationCodeLIS   4.135e-01  4.496e-01   0.920 0.357786    
## Year_fct2016:StationCodeSTTD  3.444e-01  4.849e-01   0.710 0.477618    
## Year_fct2017:StationCodeSTTD -5.748e-01  4.914e-01  -1.170 0.242151    
## Year_fct2018:StationCodeSTTD -8.724e-01  4.660e-01  -1.872 0.061252 .  
## Year_fct2019:StationCodeSTTD -1.196e+00  4.584e-01  -2.609 0.009123 ** 
## Year_fct2016:StationCodeLIB   9.423e-01  4.738e-01   1.989 0.046817 *  
## Year_fct2017:StationCodeLIB  -3.952e-01  4.620e-01  -0.855 0.392434    
## Year_fct2018:StationCodeLIB  -2.200e+00  5.051e-01  -4.356 1.36e-05 ***
## Year_fct2019:StationCodeLIB  -9.373e-01  5.127e-01  -1.828 0.067617 .  
## Year_fct2016:StationCodeRVB   3.293e-01  5.016e-01   0.657 0.511538    
## Year_fct2017:StationCodeRVB  -5.462e-01  4.829e-01  -1.131 0.258170    
## Year_fct2018:StationCodeRVB  -4.555e-01  4.624e-01  -0.985 0.324691    
## Year_fct2019:StationCodeRVB  -6.753e-01  4.697e-01  -1.438 0.150619    
## Flow:StationCodeRD22         -4.185e-04  2.542e-04  -1.646 0.099770 .  
## Flow:StationCodeI80          -1.249e-03  2.832e-04  -4.412 1.06e-05 ***
## Flow:StationCodeLIS          -2.168e-05  2.860e-04  -0.076 0.939582    
## Flow:StationCodeSTTD          1.431e-03  2.903e-04   4.930 8.61e-07 ***
## Flow:StationCodeLIB           3.575e-04  1.821e-04   1.964 0.049625 *  
## Flow:StationCodeRVB           1.822e-04  1.771e-04   1.029 0.303658    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df    F p-value    
## s(DOY) 3.675  3.675 11.3  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.819   
##   Scale est. = 0.22342   n = 3478
appraise(m_gamm_q_ar2_gam)

k.check(m_gamm_q_ar2_gam)
##        k'      edf  k-index p-value
## s(DOY)  4 3.674585 0.803638       0
draw(m_gamm_q_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_q_ar2_gam, pages = 1, all.terms = TRUE)

acf(resid_q_ar2)

Box.test(resid_q_ar2, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_q_ar2
## X-squared = 29.682, df = 20, p-value = 0.07517

The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.

# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
## 
## Parametric Terms:
##                      df      F p-value
## Year_fct              4  0.407   0.803
## Flow                  1  1.056   0.304
## StationCode           6 14.929 < 2e-16
## Year_fct:Flow         4  0.408   0.803
## Year_fct:StationCode 24  3.933 3.9e-10
## Flow:StationCode      6 16.764 < 2e-16
## 
## Approximate significance of smooth terms:
##          edf Ref.df    F p-value
## s(DOY) 3.675  3.675 11.3  <2e-16

The Year:StationCode and Flow:StationCode interactions were significant among the interaction terms of the model. Let’s remove the Year:Flow interaction term and check the model summary and diagnostics.

m_gamm_q_ar2b <- gamm(
  Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, k = 5), 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

# Pull out the residuals and the GAM model
resid_q_ar2b <- residuals(m_gamm_q_ar2b$lme, type = "normalized")
m_gamm_q_ar2b_gam <- m_gamm_q_ar2b$gam

summary(m_gamm_q_ar2b_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, 
##     k = 5)
## 
## Parametric coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   9.180e+00  2.920e-01  31.445  < 2e-16 ***
## Year_fct2016                 -3.484e-02  3.758e-01  -0.093 0.926138    
## Year_fct2017                  3.299e-01  3.584e-01   0.920 0.357500    
## Year_fct2018                  1.278e-01  3.493e-01   0.366 0.714492    
## Year_fct2019                  8.237e-02  3.503e-01   0.235 0.814093    
## StationCodeRD22               3.144e-01  3.568e-01   0.881 0.378251    
## StationCodeI80                7.079e-01  3.669e-01   1.930 0.053736 .  
## StationCodeLIS               -5.870e-01  3.494e-01  -1.680 0.093030 .  
## StationCodeSTTD              -4.804e-01  3.574e-01  -1.344 0.178970    
## StationCodeLIB               -1.239e+00  3.563e-01  -3.477 0.000514 ***
## StationCodeRVB               -1.490e+00  3.492e-01  -4.266 2.04e-05 ***
## Flow                         -1.903e-04  1.769e-04  -1.076 0.282092    
## Year_fct2016:StationCodeRD22 -4.388e-01  4.821e-01  -0.910 0.362757    
## Year_fct2017:StationCodeRD22 -5.866e-01  4.589e-01  -1.278 0.201177    
## Year_fct2018:StationCodeRD22 -4.239e-01  4.483e-01  -0.946 0.344444    
## Year_fct2019:StationCodeRD22 -1.649e-01  4.501e-01  -0.366 0.714156    
## Year_fct2016:StationCodeI80   1.310e-01  4.893e-01   0.268 0.788936    
## Year_fct2017:StationCodeI80  -3.767e-01  4.666e-01  -0.807 0.419558    
## Year_fct2018:StationCodeI80  -6.194e-01  4.562e-01  -1.358 0.174616    
## Year_fct2019:StationCodeI80   3.978e-01  4.579e-01   0.869 0.385022    
## Year_fct2016:StationCodeLIS   6.924e-01  4.697e-01   1.474 0.140536    
## Year_fct2017:StationCodeLIS   1.681e-01  4.563e-01   0.368 0.712634    
## Year_fct2018:StationCodeLIS  -1.545e-01  4.442e-01  -0.348 0.728056    
## Year_fct2019:StationCodeLIS   4.184e-01  4.493e-01   0.931 0.351866    
## Year_fct2016:StationCodeSTTD  3.480e-01  4.846e-01   0.718 0.472790    
## Year_fct2017:StationCodeSTTD -5.708e-01  4.911e-01  -1.162 0.245159    
## Year_fct2018:StationCodeSTTD -8.680e-01  4.657e-01  -1.864 0.062400 .  
## Year_fct2019:StationCodeSTTD -1.191e+00  4.581e-01  -2.600 0.009364 ** 
## Year_fct2016:StationCodeLIB   9.536e-01  4.727e-01   2.017 0.043741 *  
## Year_fct2017:StationCodeLIB  -3.846e-01  4.611e-01  -0.834 0.404321    
## Year_fct2018:StationCodeLIB  -2.191e+00  5.042e-01  -4.345 1.43e-05 ***
## Year_fct2019:StationCodeLIB  -9.096e-01  5.118e-01  -1.777 0.075641 .  
## Year_fct2016:StationCodeRVB   3.114e-01  4.819e-01   0.646 0.518156    
## Year_fct2017:StationCodeRVB  -5.575e-01  4.532e-01  -1.230 0.218786    
## Year_fct2018:StationCodeRVB  -4.597e-01  4.419e-01  -1.040 0.298248    
## Year_fct2019:StationCodeRVB  -8.461e-01  4.498e-01  -1.881 0.060059 .  
## StationCodeRD22:Flow         -4.181e-04  2.541e-04  -1.645 0.100007    
## StationCodeI80:Flow          -1.252e-03  2.831e-04  -4.424 1.00e-05 ***
## StationCodeLIS:Flow          -2.371e-05  2.859e-04  -0.083 0.933893    
## StationCodeSTTD:Flow          1.429e-03  2.902e-04   4.924 8.88e-07 ***
## StationCodeLIB:Flow           3.608e-04  1.819e-04   1.983 0.047421 *  
## StationCodeRVB:Flow           1.840e-04  1.770e-04   1.039 0.298683    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df     F p-value    
## s(DOY) 3.68   3.68 11.47  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.819   
##   Scale est. = 0.22326   n = 3478
appraise(m_gamm_q_ar2b_gam)

k.check(m_gamm_q_ar2b_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.679956 0.8026532       0
draw(m_gamm_q_ar2b_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_q_ar2b_gam, pages = 1, all.terms = TRUE)

acf(resid_q_ar2b)

Box.test(resid_q_ar2b, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_q_ar2b
## X-squared = 29.646, df = 20, p-value = 0.0758
anova(m_gamm_q_ar2b_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, 
##     k = 5)
## 
## Parametric Terms:
##                      df      F p-value
## Year_fct              4  0.406   0.805
## StationCode           6 14.873 < 2e-16
## Flow                  1  1.157   0.282
## Year_fct:StationCode 24  4.102 8.4e-11
## StationCode:Flow      6 16.865 < 2e-16
## 
## Approximate significance of smooth terms:
##         edf Ref.df     F p-value
## s(DOY) 3.68   3.68 11.47  <2e-16

This model still looks pretty good. We’ll use this one for the model selection process.

GAM Model 4: Flow and water temperature as continuous predictors

Plots

Before we add both flow and water temperature as continuous predictors in the model, let’s check if they are correlated with each other.

df_chla_wt_q %>% 
  ggplot(aes(x = Flow, y = WaterTemp)) +
  geom_point() +
  geom_smooth(formula = "y ~ x") +
  facet_wrap(vars(StationCode), scales = "free") +
  theme_bw()

Flow and water temperature appear to have no correlation.

Initial Model

We’ll try running a GAM using the model structure from the flow model and adding two-way interactions between Year:WaterTemp and Station:WaterTemp. First we’ll run the GAM without accounting for serial autocorrelation.

m_gam_q_wt <- gam(
  Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5), 
  data = df_chla_wt_q,
  method = "REML"
)

Lets look at the model summary and diagnostics:

summary(m_gam_q_wt)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * 
##     WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  10.2662032  0.2178374  47.128  < 2e-16 ***
## Year_fct2016                 -0.6484169  0.3059264  -2.120 0.034118 *  
## Year_fct2017                 -0.0609459  0.1748727  -0.349 0.727474    
## Year_fct2018                 -0.1613266  0.1680842  -0.960 0.337226    
## Year_fct2019                  0.5338094  0.1622929   3.289 0.001015 ** 
## StationCodeRD22              -1.0478545  0.1618085  -6.476 1.08e-10 ***
## StationCodeI80               -0.3198595  0.1552094  -2.061 0.039395 *  
## StationCodeLIS               -2.5319716  0.1616506 -15.663  < 2e-16 ***
## StationCodeSTTD              -1.7388040  0.1881330  -9.242  < 2e-16 ***
## StationCodeLIB               -3.9184097  0.2367461 -16.551  < 2e-16 ***
## StationCodeRVB               -3.4749646  0.1941132 -17.902  < 2e-16 ***
## Flow                         -0.0011256  0.0001014 -11.096  < 2e-16 ***
## WaterTemp                    -0.0433571  0.0092626  -4.681 2.97e-06 ***
## Year_fct2016:StationCodeRD22 -0.5347596  0.0910571  -5.873 4.69e-09 ***
## Year_fct2017:StationCodeRD22 -0.6200951  0.0889603  -6.970 3.77e-12 ***
## Year_fct2018:StationCodeRD22 -0.3517822  0.0834277  -4.217 2.54e-05 ***
## Year_fct2019:StationCodeRD22 -0.1482813  0.0833741  -1.779 0.075410 .  
## Year_fct2016:StationCodeI80   0.2012539  0.0934108   2.155 0.031270 *  
## Year_fct2017:StationCodeI80  -0.2400564  0.0912588  -2.631 0.008564 ** 
## Year_fct2018:StationCodeI80  -0.4104957  0.0855122  -4.800 1.65e-06 ***
## Year_fct2019:StationCodeI80   0.5658673  0.0856073   6.610 4.44e-11 ***
## Year_fct2016:StationCodeLIS   0.3694915  0.0892519   4.140 3.56e-05 ***
## Year_fct2017:StationCodeLIS   0.2490534  0.0878966   2.833 0.004631 ** 
## Year_fct2018:StationCodeLIS  -0.0695068  0.0825594  -0.842 0.399902    
## Year_fct2019:StationCodeLIS   0.5781116  0.0838872   6.892 6.54e-12 ***
## Year_fct2016:StationCodeSTTD  0.2714631  0.0931049   2.916 0.003572 ** 
## Year_fct2017:StationCodeSTTD -0.1673873  0.0975732  -1.716 0.086343 .  
## Year_fct2018:StationCodeSTTD -0.6929241  0.0871555  -7.950 2.50e-15 ***
## Year_fct2019:StationCodeSTTD -1.1345661  0.0858236 -13.220  < 2e-16 ***
## Year_fct2016:StationCodeLIB   0.8903568  0.0949408   9.378  < 2e-16 ***
## Year_fct2017:StationCodeLIB  -0.2315049  0.0969164  -2.389 0.016962 *  
## Year_fct2018:StationCodeLIB  -2.0080384  0.1150845 -17.448  < 2e-16 ***
## Year_fct2019:StationCodeLIB  -0.8582224  0.1196082  -7.175 8.81e-13 ***
## Year_fct2016:StationCodeRVB   0.1597391  0.1027829   1.554 0.120243    
## Year_fct2017:StationCodeRVB  -0.3753724  0.1038743  -3.614 0.000306 ***
## Year_fct2018:StationCodeRVB  -0.2920268  0.0889282  -3.284 0.001034 ** 
## Year_fct2019:StationCodeRVB  -0.6540161  0.0928199  -7.046 2.22e-12 ***
## StationCodeRD22:Flow         -0.0003430  0.0001340  -2.560 0.010519 *  
## StationCodeI80:Flow          -0.0008049  0.0001352  -5.951 2.93e-09 ***
## StationCodeLIS:Flow           0.0012604  0.0001350   9.336  < 2e-16 ***
## StationCodeSTTD:Flow          0.0030649  0.0001386  22.121  < 2e-16 ***
## StationCodeLIB:Flow           0.0013482  0.0001236  10.908  < 2e-16 ***
## StationCodeRVB:Flow           0.0011204  0.0001014  11.046  < 2e-16 ***
## Year_fct2016:WaterTemp        0.0278927  0.0123785   2.253 0.024302 *  
## Year_fct2017:WaterTemp        0.0125116  0.0069925   1.789 0.073658 .  
## Year_fct2018:WaterTemp        0.0103364  0.0068995   1.498 0.134191    
## Year_fct2019:WaterTemp       -0.0253859  0.0065403  -3.881 0.000106 ***
## StationCodeRD22:WaterTemp     0.0637778  0.0067035   9.514  < 2e-16 ***
## StationCodeI80:WaterTemp      0.0389816  0.0064179   6.074 1.39e-09 ***
## StationCodeLIS:WaterTemp      0.0857066  0.0066872  12.816  < 2e-16 ***
## StationCodeSTTD:WaterTemp     0.0496447  0.0081526   6.089 1.26e-09 ***
## StationCodeLIB:WaterTemp      0.1230615  0.0120862  10.182  < 2e-16 ***
## StationCodeRVB:WaterTemp      0.0833896  0.0086653   9.623  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df    F p-value    
## s(DOY) 3.781  3.974 53.5  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.864   Deviance explained = 86.6%
## -REML = 1615.8  Scale est. = 0.13435   n = 3478
appraise(m_gam_q_wt)

k.check(m_gam_q_wt)
##        k'     edf   k-index p-value
## s(DOY)  4 3.78136 0.9264289       0
draw(m_gam_q_wt, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gam_q_wt, pages = 1, all.terms = TRUE)

acf(residuals(m_gam_q_wt))

Box.test(residuals(m_gam_q_wt), lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuals(m_gam_q_wt)
## X-squared = 8803.9, df = 20, p-value < 2.2e-16

Model with AR() correlation structure

Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.

# Define model formula as an object
f_gam_q_wt <- as.formula(
  "Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)"
)

# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_wt <- gamm(
  f_gam_q_wt, 
  data = df_chla_wt_q,
  method = "REML"
)

m_gamm_q_wt_ar1 <- gamm(
  f_gam_q_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_wt_ar2 <- gamm(
  f_gam_q_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_wt_ar3 <- gamm(
  f_gam_q_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
  method = "REML"
)

# Compare models
anova(
  m_gamm_q_wt$lme, 
  m_gamm_q_wt_ar1$lme, 
  m_gamm_q_wt_ar2$lme, 
  m_gamm_q_wt_ar3$lme
)
##                     Model df       AIC       BIC    logLik   Test  L.Ratio
## m_gamm_q_wt$lme         1 56  3343.533  3687.292 -1615.766                
## m_gamm_q_wt_ar1$lme     2 57 -2152.733 -1802.835  1133.367 1 vs 2 5498.266
## m_gamm_q_wt_ar2$lme     3 58 -2174.391 -1818.354  1145.195 2 vs 3   23.657
## m_gamm_q_wt_ar3$lme     4 59 -2173.370 -1811.195  1145.685 3 vs 4    0.980
##                     p-value
## m_gamm_q_wt$lme            
## m_gamm_q_wt_ar1$lme  <.0001
## m_gamm_q_wt_ar2$lme  <.0001
## m_gamm_q_wt_ar3$lme  0.3222

It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.

AR(2) Model

# Pull out the residuals and the GAM model
resid_q_wt_ar2 <- residuals(m_gamm_q_wt_ar2$lme, type = "normalized")
m_gamm_q_wt_ar2_gam <- m_gamm_q_wt_ar2$gam

summary(m_gamm_q_wt_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * 
##     WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   8.894e+00  5.959e-01  14.925  < 2e-16 ***
## Year_fct2016                  1.423e+00  6.306e-01   2.256 0.024113 *  
## Year_fct2017                  4.921e-02  5.714e-01   0.086 0.931374    
## Year_fct2018                 -3.647e-01  5.778e-01  -0.631 0.527977    
## Year_fct2019                 -2.601e-02  5.490e-01  -0.047 0.962218    
## StationCodeRD22              -1.128e+00  6.141e-01  -1.838 0.066189 .  
## StationCodeI80                8.206e-01  5.957e-01   1.378 0.168447    
## StationCodeLIS               -1.605e+00  6.037e-01  -2.659 0.007879 ** 
## StationCodeSTTD              -1.252e+00  6.252e-01  -2.003 0.045228 *  
## StationCodeLIB               -9.809e-01  6.957e-01  -1.410 0.158668    
## StationCodeRVB               -2.142e+00  7.531e-01  -2.845 0.004470 ** 
## Flow                         -1.699e-04  1.767e-04  -0.961 0.336414    
## WaterTemp                     1.399e-02  2.161e-02   0.648 0.517248    
## Year_fct2016:StationCodeRD22 -6.003e-01  5.380e-01  -1.116 0.264600    
## Year_fct2017:StationCodeRD22 -6.123e-01  5.116e-01  -1.197 0.231484    
## Year_fct2018:StationCodeRD22 -3.207e-01  5.017e-01  -0.639 0.522705    
## Year_fct2019:StationCodeRD22 -1.397e-01  5.030e-01  -0.278 0.781210    
## Year_fct2016:StationCodeI80   2.047e-01  5.478e-01   0.374 0.708640    
## Year_fct2017:StationCodeI80  -2.886e-01  5.216e-01  -0.553 0.580145    
## Year_fct2018:StationCodeI80  -5.347e-01  5.115e-01  -1.045 0.295920    
## Year_fct2019:StationCodeI80   4.585e-01  5.129e-01   0.894 0.371361    
## Year_fct2016:StationCodeLIS   5.821e-01  5.244e-01   1.110 0.267053    
## Year_fct2017:StationCodeLIS   1.927e-01  5.089e-01   0.379 0.704890    
## Year_fct2018:StationCodeLIS  -1.008e-01  4.971e-01  -0.203 0.839336    
## Year_fct2019:StationCodeLIS   4.263e-01  5.023e-01   0.849 0.396085    
## Year_fct2016:StationCodeSTTD  2.585e-01  5.418e-01   0.477 0.633378    
## Year_fct2017:StationCodeSTTD -5.788e-01  5.483e-01  -1.056 0.291233    
## Year_fct2018:StationCodeSTTD -8.371e-01  5.208e-01  -1.607 0.108099    
## Year_fct2019:StationCodeSTTD -1.172e+00  5.124e-01  -2.288 0.022224 *  
## Year_fct2016:StationCodeLIB   8.689e-01  5.287e-01   1.644 0.100355    
## Year_fct2017:StationCodeLIB  -2.787e-01  5.147e-01  -0.541 0.588263    
## Year_fct2018:StationCodeLIB  -2.096e+00  5.618e-01  -3.731 0.000194 ***
## Year_fct2019:StationCodeLIB  -8.120e-01  5.685e-01  -1.428 0.153286    
## Year_fct2016:StationCodeRVB   1.962e-01  5.384e-01   0.364 0.715519    
## Year_fct2017:StationCodeRVB  -4.613e-01  5.059e-01  -0.912 0.361903    
## Year_fct2018:StationCodeRVB  -3.546e-01  4.956e-01  -0.715 0.474415    
## Year_fct2019:StationCodeRVB  -7.574e-01  5.041e-01  -1.503 0.133011    
## StationCodeRD22:Flow         -4.412e-04  2.550e-04  -1.730 0.083665 .  
## StationCodeI80:Flow          -1.205e-03  2.852e-04  -4.226 2.44e-05 ***
## StationCodeLIS:Flow          -8.722e-06  2.877e-04  -0.030 0.975817    
## StationCodeSTTD:Flow          1.426e-03  2.923e-04   4.878 1.12e-06 ***
## StationCodeLIB:Flow           3.233e-04  1.821e-04   1.775 0.075912 .  
## StationCodeRVB:Flow           1.655e-04  1.768e-04   0.936 0.349418    
## Year_fct2016:WaterTemp       -5.896e-02  2.002e-02  -2.945 0.003250 ** 
## Year_fct2017:WaterTemp        1.079e-02  1.791e-02   0.603 0.546631    
## Year_fct2018:WaterTemp        2.128e-02  1.913e-02   1.113 0.265978    
## Year_fct2019:WaterTemp        3.786e-03  1.708e-02   0.222 0.824591    
## StationCodeRD22:WaterTemp     6.651e-02  2.099e-02   3.169 0.001541 ** 
## StationCodeI80:WaterTemp     -8.905e-03  1.954e-02  -0.456 0.648553    
## StationCodeLIS:WaterTemp      4.642e-02  2.053e-02   2.261 0.023804 *  
## StationCodeSTTD:WaterTemp     3.510e-02  2.186e-02   1.606 0.108445    
## StationCodeLIB:WaterTemp     -1.635e-02  2.705e-02  -0.605 0.545485    
## StationCodeRVB:WaterTemp      2.870e-02  2.960e-02   0.970 0.332241    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.664  3.664 8.628 2.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.807   
##   Scale est. = 0.2515    n = 3478
appraise(m_gamm_q_wt_ar2_gam)

k.check(m_gamm_q_wt_ar2_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.663875 0.7855439       0
draw(m_gamm_q_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_q_wt_ar2_gam, pages = 1, all.terms = TRUE)

acf(resid_q_wt_ar2)

Box.test(resid_q_wt_ar2, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_q_wt_ar2
## X-squared = 30.519, df = 20, p-value = 0.06187

The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.

# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_wt_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * 
##     WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
## 
## Parametric Terms:
##                       df      F  p-value
## Year_fct               4  2.531 0.038608
## StationCode            6  6.652 5.20e-07
## Flow                   1  0.924 0.336414
## WaterTemp              1  0.419 0.517248
## Year_fct:StationCode  24  3.165 3.27e-07
## StationCode:Flow       6 15.268  < 2e-16
## Year_fct:WaterTemp     4  4.659 0.000945
## StationCode:WaterTemp  6  4.669 9.75e-05
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value
## s(DOY) 3.664  3.664 8.628 2.2e-06

All interaction terms are significant in this model, so we’ll use this one for the model selection process.

GAM Model 5: Flow with flow action period

Plots

Let’s look at the relationship between chlorophyll and flow for each station-flow action period combination.

df_chla_wt_q %>% 
  ggplot(aes(x = Flow, y = Chla)) +
  geom_point() +
  geom_smooth(formula = "y ~ x") +
  facet_wrap(
    vars(StationCode, FlowActionPeriod), 
    ncol = 3, 
    scales = "free", 
    labeller = labeller(.multi_line = FALSE)
  ) +
  theme_bw()

Initial Model

It’s possible that the relationship between chlorophyll and Flow differs between Flow Action Periods at each station. We’ll try to account for this by running a GAM including a two-way interaction between Year and Station, a three-way interaction between Flow, Station, and Flow Action Period, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). We won’t include the interaction between Year and Flow Action Period since it was barely significant in the categorical model and we’re primarily including the Flow Action Period to allow for different slopes for flow at each station. First we’ll run the GAM without accounting for serial autocorrelation.

m_gam_q_fa <- gam(
  Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5), 
  data = df_chla_wt_q,
  method = "REML"
)

Lets look at the model summary and diagnostics:

summary(m_gam_q_fa)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + 
##     s(DOY, k = 5)
## 
## Parametric coefficients:
##                                               Estimate Std. Error t value
## (Intercept)                                  9.0200439  0.0839467 107.450
## Year_fct2016                                -0.0752673  0.0817958  -0.920
## Year_fct2017                                 0.1772134  0.0707049   2.506
## Year_fct2018                                -0.0051101  0.0673845  -0.076
## Year_fct2019                                -0.0504808  0.0701749  -0.719
## StationCodeRD22                              0.6048920  0.1021743   5.920
## StationCodeI80                               0.7075466  0.1053335   6.717
## StationCodeLIS                              -0.2534964  0.1008109  -2.515
## StationCodeSTTD                             -0.4888704  0.1077499  -4.537
## StationCodeLIB                              -0.5233575  0.1654118  -3.164
## StationCodeRVB                              -1.3296625  0.1140411 -11.660
## Flow                                         0.0003845  0.0006336   0.607
## FlowActionPeriodDuring                       0.2206777  0.0876376   2.518
## FlowActionPeriodAfter                        0.8225799  0.0753879  10.911
## Year_fct2016:StationCodeRD22                -0.3499073  0.0955876  -3.661
## Year_fct2017:StationCodeRD22                -0.5693333  0.0856991  -6.643
## Year_fct2018:StationCodeRD22                -0.3195979  0.0814168  -3.925
## Year_fct2019:StationCodeRD22                -0.0887040  0.0850886  -1.042
## Year_fct2016:StationCodeI80                  0.3049614  0.0984091   3.099
## Year_fct2017:StationCodeI80                 -0.2386526  0.0866819  -2.753
## Year_fct2018:StationCodeI80                 -0.3895701  0.0840942  -4.633
## Year_fct2019:StationCodeI80                  0.6211290  0.0862887   7.198
## Year_fct2016:StationCodeLIS                  0.3767691  0.0932856   4.039
## Year_fct2017:StationCodeLIS                  0.2509619  0.0859328   2.920
## Year_fct2018:StationCodeLIS                  0.0062195  0.0810131   0.077
## Year_fct2019:StationCodeLIS                  0.6367041  0.0862458   7.382
## Year_fct2016:StationCodeSTTD                 0.2375629  0.0960514   2.473
## Year_fct2017:StationCodeSTTD                -0.0453964  0.0976543  -0.465
## Year_fct2018:StationCodeSTTD                -0.5512862  0.0857425  -6.430
## Year_fct2019:StationCodeSTTD                -1.0150911  0.0885178 -11.468
## Year_fct2016:StationCodeLIB                  0.8861875  0.0923868   9.592
## Year_fct2017:StationCodeLIB                 -0.2709091  0.0930803  -2.910
## Year_fct2018:StationCodeLIB                 -2.0402939  0.1103463 -18.490
## Year_fct2019:StationCodeLIB                 -1.0231592  0.1184799  -8.636
## Year_fct2016:StationCodeRVB                  0.1576659  0.1008718   1.563
## Year_fct2017:StationCodeRVB                 -0.4120981  0.0981257  -4.200
## Year_fct2018:StationCodeRVB                 -0.2930854  0.0862916  -3.396
## Year_fct2019:StationCodeRVB                 -0.6454502  0.0914059  -7.061
## StationCodeRD22:Flow                         0.0003861  0.0014058   0.275
## StationCodeI80:Flow                          0.0018616  0.0018957   0.982
## StationCodeLIS:Flow                          0.0030944  0.0009665   3.202
## StationCodeSTTD:Flow                         0.0023417  0.0011000   2.129
## StationCodeLIB:Flow                          0.0001091  0.0006393   0.171
## StationCodeRVB:Flow                         -0.0003875  0.0006337  -0.612
## Flow:FlowActionPeriodDuring                 -0.0013504  0.0006478  -2.084
## Flow:FlowActionPeriodAfter                  -0.0083453  0.0011337  -7.361
## StationCodeRD22:FlowActionPeriodDuring      -0.2621276  0.1125150  -2.330
## StationCodeI80:FlowActionPeriodDuring       -0.1046615  0.1170017  -0.895
## StationCodeLIS:FlowActionPeriodDuring        0.3486851  0.1118607   3.117
## StationCodeSTTD:FlowActionPeriodDuring       0.4101368  0.1221588   3.357
## StationCodeLIB:FlowActionPeriodDuring       -0.3116617  0.1631523  -1.910
## StationCodeRVB:FlowActionPeriodDuring       -0.2557270  0.1403998  -1.821
## StationCodeRD22:FlowActionPeriodAfter       -0.6738249  0.0876828  -7.685
## StationCodeI80:FlowActionPeriodAfter        -0.4696770  0.0926079  -5.072
## StationCodeLIS:FlowActionPeriodAfter        -0.6118818  0.0862036  -7.098
## StationCodeSTTD:FlowActionPeriodAfter       -0.4414052  0.0972779  -4.538
## StationCodeLIB:FlowActionPeriodAfter        -1.9012296  0.1443312 -13.173
## StationCodeRVB:FlowActionPeriodAfter        -0.9357160  0.1208432  -7.743
## StationCodeRD22:Flow:FlowActionPeriodDuring -0.0007004  0.0014165  -0.494
## StationCodeI80:Flow:FlowActionPeriodDuring  -0.0028106  0.0019000  -1.479
## StationCodeLIS:Flow:FlowActionPeriodDuring  -0.0035976  0.0009869  -3.645
## StationCodeSTTD:Flow:FlowActionPeriodDuring -0.0007149  0.0011208  -0.638
## StationCodeLIB:Flow:FlowActionPeriodDuring   0.0011745  0.0006559   1.791
## StationCodeRVB:Flow:FlowActionPeriodDuring   0.0013536  0.0006479   2.089
## StationCodeRD22:Flow:FlowActionPeriodAfter  -0.0017581  0.0023433  -0.750
## StationCodeI80:Flow:FlowActionPeriodAfter   -0.0055463  0.0025798  -2.150
## StationCodeLIS:Flow:FlowActionPeriodAfter    0.0125081  0.0015103   8.282
## StationCodeSTTD:Flow:FlowActionPeriodAfter   0.0118661  0.0016282   7.288
## StationCodeLIB:Flow:FlowActionPeriodAfter    0.0076975  0.0011391   6.758
## StationCodeRVB:Flow:FlowActionPeriodAfter    0.0083444  0.0011336   7.361
##                                             Pr(>|t|)    
## (Intercept)                                  < 2e-16 ***
## Year_fct2016                                0.357541    
## Year_fct2017                                0.012244 *  
## Year_fct2018                                0.939555    
## Year_fct2019                                0.471970    
## StationCodeRD22                             3.53e-09 ***
## StationCodeI80                              2.16e-11 ***
## StationCodeLIS                              0.011963 *  
## StationCodeSTTD                             5.90e-06 ***
## StationCodeLIB                              0.001570 ** 
## StationCodeRVB                               < 2e-16 ***
## Flow                                        0.543936    
## FlowActionPeriodDuring                      0.011846 *  
## FlowActionPeriodAfter                        < 2e-16 ***
## Year_fct2016:StationCodeRD22                0.000255 ***
## Year_fct2017:StationCodeRD22                3.55e-11 ***
## Year_fct2018:StationCodeRD22                8.83e-05 ***
## Year_fct2019:StationCodeRD22                0.297259    
## Year_fct2016:StationCodeI80                 0.001958 ** 
## Year_fct2017:StationCodeI80                 0.005933 ** 
## Year_fct2018:StationCodeI80                 3.75e-06 ***
## Year_fct2019:StationCodeI80                 7.47e-13 ***
## Year_fct2016:StationCodeLIS                 5.49e-05 ***
## Year_fct2017:StationCodeLIS                 0.003518 ** 
## Year_fct2018:StationCodeLIS                 0.938810    
## Year_fct2019:StationCodeLIS                 1.94e-13 ***
## Year_fct2016:StationCodeSTTD                0.013436 *  
## Year_fct2017:StationCodeSTTD                0.642055    
## Year_fct2018:StationCodeSTTD                1.46e-10 ***
## Year_fct2019:StationCodeSTTD                 < 2e-16 ***
## Year_fct2016:StationCodeLIB                  < 2e-16 ***
## Year_fct2017:StationCodeLIB                 0.003632 ** 
## Year_fct2018:StationCodeLIB                  < 2e-16 ***
## Year_fct2019:StationCodeLIB                  < 2e-16 ***
## Year_fct2016:StationCodeRVB                 0.118138    
## Year_fct2017:StationCodeRVB                 2.74e-05 ***
## Year_fct2018:StationCodeRVB                 0.000690 ***
## Year_fct2019:StationCodeRVB                 1.99e-12 ***
## StationCodeRD22:Flow                        0.783602    
## StationCodeI80:Flow                         0.326154    
## StationCodeLIS:Flow                         0.001378 ** 
## StationCodeSTTD:Flow                        0.033342 *  
## StationCodeLIB:Flow                         0.864475    
## StationCodeRVB:Flow                         0.540909    
## Flow:FlowActionPeriodDuring                 0.037198 *  
## Flow:FlowActionPeriodAfter                  2.27e-13 ***
## StationCodeRD22:FlowActionPeriodDuring      0.019880 *  
## StationCodeI80:FlowActionPeriodDuring       0.371102    
## StationCodeLIS:FlowActionPeriodDuring       0.001841 ** 
## StationCodeSTTD:FlowActionPeriodDuring      0.000795 ***
## StationCodeLIB:FlowActionPeriodDuring       0.056185 .  
## StationCodeRVB:FlowActionPeriodDuring       0.068631 .  
## StationCodeRD22:FlowActionPeriodAfter       1.99e-14 ***
## StationCodeI80:FlowActionPeriodAfter        4.15e-07 ***
## StationCodeLIS:FlowActionPeriodAfter        1.53e-12 ***
## StationCodeSTTD:FlowActionPeriodAfter       5.89e-06 ***
## StationCodeLIB:FlowActionPeriodAfter         < 2e-16 ***
## StationCodeRVB:FlowActionPeriodAfter        1.27e-14 ***
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.621005    
## StationCodeI80:Flow:FlowActionPeriodDuring  0.139151    
## StationCodeLIS:Flow:FlowActionPeriodDuring  0.000271 ***
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.523617    
## StationCodeLIB:Flow:FlowActionPeriodDuring  0.073440 .  
## StationCodeRVB:Flow:FlowActionPeriodDuring  0.036769 *  
## StationCodeRD22:Flow:FlowActionPeriodAfter  0.453163    
## StationCodeI80:Flow:FlowActionPeriodAfter   0.031634 *  
## StationCodeLIS:Flow:FlowActionPeriodAfter    < 2e-16 ***
## StationCodeSTTD:Flow:FlowActionPeriodAfter  3.90e-13 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter   1.64e-11 ***
## StationCodeRVB:Flow:FlowActionPeriodAfter   2.28e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.863  3.989 47.56  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.88   Deviance explained = 88.2%
## -REML = 1476.1  Scale est. = 0.11868   n = 3478
appraise(m_gam_q_fa)

k.check(m_gam_q_fa)
##        k'      edf   k-index p-value
## s(DOY)  4 3.862781 0.9687532  0.0375
draw(m_gam_q_fa, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gam_q_fa, pages = 1, all.terms = TRUE)

acf(residuals(m_gam_q_fa))

Box.test(residuals(m_gam_q_fa), lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuals(m_gam_q_fa)
## X-squared = 7116.1, df = 20, p-value < 2.2e-16

Model with AR() correlation structure

Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.

# Define model formula as an object
f_gam_q_fa <- as.formula("Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5)")

# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_fa <- gamm(
  f_gam_q_fa, 
  data = df_chla_wt_q,
  method = "REML"
)

m_gamm_q_fa_ar1 <- gamm(
  f_gam_q_fa, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_fa_ar2 <- gamm(
  f_gam_q_fa, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_fa_ar3 <- gamm(
  f_gam_q_fa, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
  method = "REML"
)

# Compare models
anova(
  m_gamm_q_fa$lme, 
  m_gamm_q_fa_ar1$lme, 
  m_gamm_q_fa_ar2$lme, 
  m_gamm_q_fa_ar3$lme
)
##                     Model df       AIC       BIC    logLik   Test  L.Ratio
## m_gamm_q_fa$lme         1 73  3098.134  3545.886 -1476.067                
## m_gamm_q_fa_ar1$lme     2 74 -1985.717 -1531.832  1066.859 1 vs 2 5085.851
## m_gamm_q_fa_ar2$lme     3 75 -2005.949 -1545.930  1077.975 2 vs 3   22.232
## m_gamm_q_fa_ar3$lme     4 76 -2004.214 -1538.061  1078.107 3 vs 4    0.265
##                     p-value
## m_gamm_q_fa$lme            
## m_gamm_q_fa_ar1$lme  <.0001
## m_gamm_q_fa_ar2$lme  <.0001
## m_gamm_q_fa_ar3$lme  0.6067

Again, it looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.

AR(2) Model

# Pull out the residuals and the GAM model
resid_q_fa_ar2 <- residuals(m_gamm_q_fa_ar2$lme, type = "normalized")
m_gamm_q_fa_ar2_gam <- m_gamm_q_fa_ar2$gam

summary(m_gamm_q_fa_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + 
##     s(DOY, k = 5)
## 
## Parametric coefficients:
##                                               Estimate Std. Error t value
## (Intercept)                                  9.140e+00  3.300e-01  27.694
## Year_fct2016                                -8.591e-02  4.210e-01  -0.204
## Year_fct2017                                 3.071e-01  3.994e-01   0.769
## Year_fct2018                                 1.029e-01  3.897e-01   0.264
## Year_fct2019                                 6.427e-02  3.912e-01   0.164
## StationCodeRD22                              2.084e-01  4.107e-01   0.508
## StationCodeI80                               3.448e-01  4.231e-01   0.815
## StationCodeLIS                              -6.571e-01  3.993e-01  -1.646
## StationCodeSTTD                             -2.267e-01  4.099e-01  -0.553
## StationCodeLIB                              -9.494e-01  4.148e-01  -2.289
## StationCodeRVB                              -1.386e+00  4.029e-01  -3.441
## Flow                                        -1.686e-04  3.656e-04  -0.461
## FlowActionPeriodDuring                       8.981e-02  1.096e-01   0.819
## FlowActionPeriodAfter                        1.447e-01  1.313e-01   1.102
## Year_fct2016:StationCodeRD22                -4.345e-01  5.402e-01  -0.804
## Year_fct2017:StationCodeRD22                -5.513e-01  5.127e-01  -1.075
## Year_fct2018:StationCodeRD22                -3.939e-01  5.015e-01  -0.785
## Year_fct2019:StationCodeRD22                -1.286e-01  5.039e-01  -0.255
## Year_fct2016:StationCodeI80                  1.236e-01  5.486e-01   0.225
## Year_fct2017:StationCodeI80                 -3.138e-01  5.212e-01  -0.602
## Year_fct2018:StationCodeI80                 -5.976e-01  5.102e-01  -1.171
## Year_fct2019:StationCodeI80                  4.697e-01  5.125e-01   0.916
## Year_fct2016:StationCodeLIS                  7.761e-01  5.260e-01   1.475
## Year_fct2017:StationCodeLIS                  2.046e-01  5.104e-01   0.401
## Year_fct2018:StationCodeLIS                 -1.434e-01  4.971e-01  -0.288
## Year_fct2019:StationCodeLIS                  4.044e-01  5.029e-01   0.804
## Year_fct2016:StationCodeSTTD                 3.971e-01  5.426e-01   0.732
## Year_fct2017:StationCodeSTTD                -6.679e-01  5.490e-01  -1.217
## Year_fct2018:StationCodeSTTD                -9.501e-01  5.208e-01  -1.824
## Year_fct2019:StationCodeSTTD                -1.144e+00  5.126e-01  -2.232
## Year_fct2016:StationCodeLIB                  9.970e-01  5.295e-01   1.883
## Year_fct2017:StationCodeLIB                 -4.086e-01  5.154e-01  -0.793
## Year_fct2018:StationCodeLIB                 -2.176e+00  5.607e-01  -3.881
## Year_fct2019:StationCodeLIB                 -9.519e-01  5.694e-01  -1.672
## Year_fct2016:StationCodeRVB                  3.678e-01  5.390e-01   0.682
## Year_fct2017:StationCodeRVB                 -5.412e-01  5.064e-01  -1.069
## Year_fct2018:StationCodeRVB                 -4.418e-01  4.948e-01  -0.893
## Year_fct2019:StationCodeRVB                 -8.452e-01  5.035e-01  -1.678
## StationCodeRD22:Flow                         6.557e-04  1.189e-03   0.552
## StationCodeI80:Flow                          4.714e-03  1.602e-03   2.943
## StationCodeLIS:Flow                         -8.155e-04  8.945e-04  -0.912
## StationCodeSTTD:Flow                        -1.186e-04  1.026e-03  -0.116
## StationCodeLIB:Flow                          4.674e-04  3.725e-04   1.255
## StationCodeRVB:Flow                          1.588e-04  3.657e-04   0.434
## Flow:FlowActionPeriodDuring                 -5.381e-05  3.758e-04  -0.143
## Flow:FlowActionPeriodAfter                   7.104e-05  7.263e-04   0.098
## StationCodeRD22:FlowActionPeriodDuring       3.777e-02  1.509e-01   0.250
## StationCodeI80:FlowActionPeriodDuring        3.444e-01  1.540e-01   2.236
## StationCodeLIS:FlowActionPeriodDuring       -1.107e-01  1.328e-01  -0.834
## StationCodeSTTD:FlowActionPeriodDuring       7.934e-03  1.340e-01   0.059
## StationCodeLIB:FlowActionPeriodDuring       -4.421e-01  1.783e-01  -2.480
## StationCodeRVB:FlowActionPeriodDuring       -1.782e-01  1.846e-01  -0.965
## StationCodeRD22:FlowActionPeriodAfter        8.258e-02  1.954e-01   0.423
## StationCodeI80:FlowActionPeriodAfter         4.533e-01  1.984e-01   2.285
## StationCodeLIS:FlowActionPeriodAfter         1.976e-02  1.635e-01   0.121
## StationCodeSTTD:FlowActionPeriodAfter       -4.926e-01  1.658e-01  -2.972
## StationCodeLIB:FlowActionPeriodAfter        -4.578e-01  1.974e-01  -2.320
## StationCodeRVB:FlowActionPeriodAfter        -2.304e-01  1.999e-01  -1.153
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.040e-03  1.172e-03  -0.887
## StationCodeI80:Flow:FlowActionPeriodDuring  -5.914e-03  1.569e-03  -3.771
## StationCodeLIS:Flow:FlowActionPeriodDuring   1.118e-03  9.342e-04   1.197
## StationCodeSTTD:Flow:FlowActionPeriodDuring  9.933e-04  1.061e-03   0.936
## StationCodeLIB:Flow:FlowActionPeriodDuring  -1.879e-04  3.889e-04  -0.483
## StationCodeRVB:Flow:FlowActionPeriodDuring   6.217e-05  3.761e-04   0.165
## StationCodeRD22:Flow:FlowActionPeriodAfter  -4.517e-05  2.549e-03  -0.018
## StationCodeI80:Flow:FlowActionPeriodAfter   -8.338e-03  2.762e-03  -3.019
## StationCodeLIS:Flow:FlowActionPeriodAfter   -4.819e-05  1.269e-03  -0.038
## StationCodeSTTD:Flow:FlowActionPeriodAfter   4.786e-03  1.381e-03   3.465
## StationCodeLIB:Flow:FlowActionPeriodAfter   -2.209e-04  7.319e-04  -0.302
## StationCodeRVB:Flow:FlowActionPeriodAfter   -6.675e-05  7.263e-04  -0.092
##                                             Pr(>|t|)    
## (Intercept)                                  < 2e-16 ***
## Year_fct2016                                0.838299    
## Year_fct2017                                0.441995    
## Year_fct2018                                0.791690    
## Year_fct2019                                0.869520    
## StationCodeRD22                             0.611816    
## StationCodeI80                              0.415252    
## StationCodeLIS                              0.099916 .  
## StationCodeSTTD                             0.580352    
## StationCodeLIB                              0.022136 *  
## StationCodeRVB                              0.000587 ***
## Flow                                        0.644739    
## FlowActionPeriodDuring                      0.412618    
## FlowActionPeriodAfter                       0.270588    
## Year_fct2016:StationCodeRD22                0.421258    
## Year_fct2017:StationCodeRD22                0.282374    
## Year_fct2018:StationCodeRD22                0.432237    
## Year_fct2019:StationCodeRD22                0.798634    
## Year_fct2016:StationCodeI80                 0.821794    
## Year_fct2017:StationCodeI80                 0.547132    
## Year_fct2018:StationCodeI80                 0.241558    
## Year_fct2019:StationCodeI80                 0.359519    
## Year_fct2016:StationCodeLIS                 0.140224    
## Year_fct2017:StationCodeLIS                 0.688545    
## Year_fct2018:StationCodeLIS                 0.772996    
## Year_fct2019:StationCodeLIS                 0.421423    
## Year_fct2016:StationCodeSTTD                0.464298    
## Year_fct2017:StationCodeSTTD                0.223845    
## Year_fct2018:StationCodeSTTD                0.068178 .  
## Year_fct2019:StationCodeSTTD                0.025692 *  
## Year_fct2016:StationCodeLIB                 0.059765 .  
## Year_fct2017:StationCodeLIB                 0.427980    
## Year_fct2018:StationCodeLIB                 0.000106 ***
## Year_fct2019:StationCodeLIB                 0.094640 .  
## Year_fct2016:StationCodeRVB                 0.495041    
## Year_fct2017:StationCodeRVB                 0.285267    
## Year_fct2018:StationCodeRVB                 0.371985    
## Year_fct2019:StationCodeRVB                 0.093350 .  
## StationCodeRD22:Flow                        0.581215    
## StationCodeI80:Flow                         0.003268 ** 
## StationCodeLIS:Flow                         0.362050    
## StationCodeSTTD:Flow                        0.908046    
## StationCodeLIB:Flow                         0.209715    
## StationCodeRVB:Flow                         0.664224    
## Flow:FlowActionPeriodDuring                 0.886142    
## Flow:FlowActionPeriodAfter                  0.922079    
## StationCodeRD22:FlowActionPeriodDuring      0.802365    
## StationCodeI80:FlowActionPeriodDuring       0.025409 *  
## StationCodeLIS:FlowActionPeriodDuring       0.404406    
## StationCodeSTTD:FlowActionPeriodDuring      0.952780    
## StationCodeLIB:FlowActionPeriodDuring       0.013178 *  
## StationCodeRVB:FlowActionPeriodDuring       0.334379    
## StationCodeRD22:FlowActionPeriodAfter       0.672656    
## StationCodeI80:FlowActionPeriodAfter        0.022380 *  
## StationCodeLIS:FlowActionPeriodAfter        0.903816    
## StationCodeSTTD:FlowActionPeriodAfter       0.002984 ** 
## StationCodeLIB:FlowActionPeriodAfter        0.020421 *  
## StationCodeRVB:FlowActionPeriodAfter        0.249071    
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.374929    
## StationCodeI80:Flow:FlowActionPeriodDuring  0.000166 ***
## StationCodeLIS:Flow:FlowActionPeriodDuring  0.231434    
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.349200    
## StationCodeLIB:Flow:FlowActionPeriodDuring  0.629007    
## StationCodeRVB:Flow:FlowActionPeriodDuring  0.868712    
## StationCodeRD22:Flow:FlowActionPeriodAfter  0.985860    
## StationCodeI80:Flow:FlowActionPeriodAfter   0.002555 ** 
## StationCodeLIS:Flow:FlowActionPeriodAfter   0.969721    
## StationCodeSTTD:Flow:FlowActionPeriodAfter  0.000537 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter   0.762798    
## StationCodeRVB:Flow:FlowActionPeriodAfter   0.926785    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value    
## s(DOY) 3.651  3.651 9.374 1.24e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.806   
##   Scale est. = 0.25102   n = 3478
appraise(m_gamm_q_fa_ar2_gam)

k.check(m_gamm_q_fa_ar2_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.650986 0.8149469       0
draw(m_gamm_q_fa_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_q_fa_ar2_gam, pages = 1, all.terms = TRUE)

acf(resid_q_fa_ar2)

Box.test(resid_q_fa_ar2, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_q_fa_ar2
## X-squared = 30.35, df = 20, p-value = 0.06438

The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.

# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_fa_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + 
##     s(DOY, k = 5)
## 
## Parametric Terms:
##                                   df     F  p-value
## Year_fct                           4 0.342 0.849425
## StationCode                        6 6.770 3.80e-07
## Flow                               1 0.213 0.644739
## FlowActionPeriod                   2 0.607 0.544928
## Year_fct:StationCode              24 3.414 3.87e-08
## StationCode:Flow                   6 4.543 0.000135
## Flow:FlowActionPeriod              2 0.026 0.974524
## StationCode:FlowActionPeriod      12 5.610 1.29e-09
## StationCode:Flow:FlowActionPeriod 12 5.040 2.22e-08
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value
## s(DOY) 3.651  3.651 9.374 1.24e-06

Despite its complexity, we’ll use this model for the model selection process.

GAM Model 6: Flow with flow action period and water temperature

Initial Model

Finally, we’ll try running a GAM using the model structure from the flow with flow action period model and adding two-way interactions between Year:WaterTemp and Station:WaterTemp. First we’ll run the GAM without accounting for serial autocorrelation.

m_gam_q_fa_wt <- gam(
  Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5), 
  data = df_chla_wt_q,
  method = "REML"
)

Lets look at the model summary and diagnostics:

summary(m_gam_q_fa_wt)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + 
##     Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                               Estimate Std. Error t value
## (Intercept)                                  9.547e+00  2.926e-01  32.628
## Year_fct2016                                -3.523e-01  3.056e-01  -1.153
## Year_fct2017                                -5.764e-03  1.687e-01  -0.034
## Year_fct2018                                -3.353e-01  1.608e-01  -2.086
## Year_fct2019                                 4.571e-01  1.558e-01   2.933
## StationCodeRD22                             -4.980e-01  3.355e-01  -1.484
## StationCodeI80                               7.357e-01  3.224e-01   2.281
## StationCodeLIS                              -2.400e+00  2.969e-01  -8.084
## StationCodeSTTD                             -1.276e+00  3.227e-01  -3.955
## StationCodeLIB                              -1.247e+00  4.001e-01  -3.117
## StationCodeRVB                              -2.148e+00  3.514e-01  -6.113
## Flow                                        -7.271e-04  6.433e-04  -1.130
## FlowActionPeriodDuring                       3.434e-02  9.046e-02   0.380
## FlowActionPeriodAfter                        4.980e-01  1.057e-01   4.710
## WaterTemp                                   -1.568e-02  1.146e-02  -1.368
## Year_fct2016:StationCodeRD22                -4.571e-01  1.005e-01  -4.548
## Year_fct2017:StationCodeRD22                -5.739e-01  8.460e-02  -6.783
## Year_fct2018:StationCodeRD22                -2.862e-01  8.045e-02  -3.558
## Year_fct2019:StationCodeRD22                -6.503e-02  8.367e-02  -0.777
## Year_fct2016:StationCodeI80                  3.141e-01  1.040e-01   3.022
## Year_fct2017:StationCodeI80                 -1.949e-01  8.623e-02  -2.260
## Year_fct2018:StationCodeI80                 -3.712e-01  8.364e-02  -4.438
## Year_fct2019:StationCodeI80                  6.428e-01  8.528e-02   7.538
## Year_fct2016:StationCodeLIS                  1.408e-01  9.922e-02   1.419
## Year_fct2017:StationCodeLIS                  1.859e-01  8.511e-02   2.185
## Year_fct2018:StationCodeLIS                  1.916e-02  7.961e-02   0.241
## Year_fct2019:StationCodeLIS                  6.092e-01  8.473e-02   7.190
## Year_fct2016:StationCodeSTTD                 1.687e-01  1.025e-01   1.647
## Year_fct2017:StationCodeSTTD                -5.588e-02  9.667e-02  -0.578
## Year_fct2018:StationCodeSTTD                -5.559e-01  8.433e-02  -6.592
## Year_fct2019:StationCodeSTTD                -1.018e+00  8.731e-02 -11.662
## Year_fct2016:StationCodeLIB                  8.683e-01  1.027e-01   8.458
## Year_fct2017:StationCodeLIB                 -2.345e-01  9.240e-02  -2.538
## Year_fct2018:StationCodeLIB                 -1.994e+00  1.094e-01 -18.217
## Year_fct2019:StationCodeLIB                 -9.703e-01  1.172e-01  -8.279
## Year_fct2016:StationCodeRVB                  1.142e-01  1.092e-01   1.046
## Year_fct2017:StationCodeRVB                 -3.724e-01  9.999e-02  -3.724
## Year_fct2018:StationCodeRVB                 -2.509e-01  8.781e-02  -2.857
## Year_fct2019:StationCodeRVB                 -6.336e-01  9.285e-02  -6.824
## StationCodeRD22:Flow                         6.978e-04  1.383e-03   0.505
## StationCodeI80:Flow                          4.347e-04  1.870e-03   0.232
## StationCodeLIS:Flow                          5.075e-03  9.885e-04   5.134
## StationCodeSTTD:Flow                         3.462e-03  1.099e-03   3.149
## StationCodeLIB:Flow                          1.214e-03  6.495e-04   1.870
## StationCodeRVB:Flow                          7.217e-04  6.433e-04   1.122
## Flow:FlowActionPeriodDuring                 -1.792e-04  6.558e-04  -0.273
## Flow:FlowActionPeriodAfter                  -6.344e-03  1.178e-03  -5.385
## StationCodeRD22:FlowActionPeriodDuring      -1.234e-01  1.156e-01  -1.068
## StationCodeI80:FlowActionPeriodDuring       -8.921e-02  1.197e-01  -0.746
## StationCodeLIS:FlowActionPeriodDuring        5.136e-01  1.132e-01   4.539
## StationCodeSTTD:FlowActionPeriodDuring       5.521e-01  1.237e-01   4.463
## StationCodeLIB:FlowActionPeriodDuring       -1.717e-01  1.636e-01  -1.050
## StationCodeRVB:FlowActionPeriodDuring       -1.081e-01  1.414e-01  -0.764
## StationCodeRD22:FlowActionPeriodAfter       -2.807e-01  1.429e-01  -1.965
## StationCodeI80:FlowActionPeriodAfter        -5.002e-01  1.462e-01  -3.422
## StationCodeLIS:FlowActionPeriodAfter        -4.052e-02  1.220e-01  -0.332
## StationCodeSTTD:FlowActionPeriodAfter       -1.510e-01  1.302e-01  -1.160
## StationCodeLIB:FlowActionPeriodAfter        -1.600e+00  1.802e-01  -8.878
## StationCodeRVB:FlowActionPeriodAfter        -6.906e-01  1.593e-01  -4.336
## Year_fct2016:WaterTemp                       1.554e-02  1.201e-02   1.294
## Year_fct2017:WaterTemp                       7.768e-03  6.761e-03   1.149
## Year_fct2018:WaterTemp                       1.534e-02  6.559e-03   2.339
## Year_fct2019:WaterTemp                      -2.449e-02  6.218e-03  -3.938
## StationCodeRD22:WaterTemp                    4.178e-02  1.249e-02   3.345
## StationCodeI80:WaterTemp                    -2.552e-03  1.187e-02  -0.215
## StationCodeLIS:WaterTemp                     8.755e-02  1.119e-02   7.823
## StationCodeSTTD:WaterTemp                    2.874e-02  1.227e-02   2.343
## StationCodeLIB:WaterTemp                     2.504e-02  1.611e-02   1.555
## StationCodeRVB:WaterTemp                     3.059e-02  1.375e-02   2.224
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.012e-03  1.393e-03  -0.727
## StationCodeI80:Flow:FlowActionPeriodDuring  -1.383e-03  1.874e-03  -0.738
## StationCodeLIS:Flow:FlowActionPeriodDuring  -5.526e-03  1.007e-03  -5.491
## StationCodeSTTD:Flow:FlowActionPeriodDuring -1.822e-03  1.117e-03  -1.631
## StationCodeLIB:Flow:FlowActionPeriodDuring  -1.014e-05  6.641e-04  -0.015
## StationCodeRVB:Flow:FlowActionPeriodDuring   1.844e-04  6.559e-04   0.281
## StationCodeRD22:Flow:FlowActionPeriodAfter  -4.200e-03  2.534e-03  -1.657
## StationCodeI80:Flow:FlowActionPeriodAfter   -1.946e-03  2.783e-03  -0.699
## StationCodeLIS:Flow:FlowActionPeriodAfter    8.734e-03  1.560e-03   5.598
## StationCodeSTTD:Flow:FlowActionPeriodAfter   1.013e-02  1.662e-03   6.094
## StationCodeLIB:Flow:FlowActionPeriodAfter    5.708e-03  1.184e-03   4.820
## StationCodeRVB:Flow:FlowActionPeriodAfter    6.351e-03  1.178e-03   5.389
##                                             Pr(>|t|)    
## (Intercept)                                  < 2e-16 ***
## Year_fct2016                                0.249091    
## Year_fct2017                                0.972744    
## Year_fct2018                                0.037076 *  
## Year_fct2019                                0.003375 ** 
## StationCodeRD22                             0.137856    
## StationCodeI80                              0.022584 *  
## StationCodeLIS                              8.60e-16 ***
## StationCodeSTTD                             7.82e-05 ***
## StationCodeLIB                              0.001842 ** 
## StationCodeRVB                              1.09e-09 ***
## Flow                                        0.258469    
## FlowActionPeriodDuring                      0.704248    
## FlowActionPeriodAfter                       2.58e-06 ***
## WaterTemp                                   0.171506    
## Year_fct2016:StationCodeRD22                5.60e-06 ***
## Year_fct2017:StationCodeRD22                1.38e-11 ***
## Year_fct2018:StationCodeRD22                0.000379 ***
## Year_fct2019:StationCodeRD22                0.437081    
## Year_fct2016:StationCodeI80                 0.002533 ** 
## Year_fct2017:StationCodeI80                 0.023879 *  
## Year_fct2018:StationCodeI80                 9.35e-06 ***
## Year_fct2019:StationCodeI80                 6.10e-14 ***
## Year_fct2016:StationCodeLIS                 0.155939    
## Year_fct2017:StationCodeLIS                 0.028987 *  
## Year_fct2018:StationCodeLIS                 0.809791    
## Year_fct2019:StationCodeLIS                 7.93e-13 ***
## Year_fct2016:StationCodeSTTD                0.099692 .  
## Year_fct2017:StationCodeSTTD                0.563289    
## Year_fct2018:StationCodeSTTD                5.03e-11 ***
## Year_fct2019:StationCodeSTTD                 < 2e-16 ***
## Year_fct2016:StationCodeLIB                  < 2e-16 ***
## Year_fct2017:StationCodeLIB                 0.011181 *  
## Year_fct2018:StationCodeLIB                  < 2e-16 ***
## Year_fct2019:StationCodeLIB                  < 2e-16 ***
## Year_fct2016:StationCodeRVB                 0.295425    
## Year_fct2017:StationCodeRVB                 0.000199 ***
## Year_fct2018:StationCodeRVB                 0.004298 ** 
## Year_fct2019:StationCodeRVB                 1.05e-11 ***
## StationCodeRD22:Flow                        0.613905    
## StationCodeI80:Flow                         0.816208    
## StationCodeLIS:Flow                         2.99e-07 ***
## StationCodeSTTD:Flow                        0.001652 ** 
## StationCodeLIB:Flow                         0.061637 .  
## StationCodeRVB:Flow                         0.261987    
## Flow:FlowActionPeriodDuring                 0.784664    
## Flow:FlowActionPeriodAfter                  7.73e-08 ***
## StationCodeRD22:FlowActionPeriodDuring      0.285703    
## StationCodeI80:FlowActionPeriodDuring       0.455987    
## StationCodeLIS:FlowActionPeriodDuring       5.86e-06 ***
## StationCodeSTTD:FlowActionPeriodDuring      8.33e-06 ***
## StationCodeLIB:FlowActionPeriodDuring       0.293788    
## StationCodeRVB:FlowActionPeriodDuring       0.444642    
## StationCodeRD22:FlowActionPeriodAfter       0.049498 *  
## StationCodeI80:FlowActionPeriodAfter        0.000629 ***
## StationCodeLIS:FlowActionPeriodAfter        0.739835    
## StationCodeSTTD:FlowActionPeriodAfter       0.246237    
## StationCodeLIB:FlowActionPeriodAfter         < 2e-16 ***
## StationCodeRVB:FlowActionPeriodAfter        1.49e-05 ***
## Year_fct2016:WaterTemp                      0.195712    
## Year_fct2017:WaterTemp                      0.250655    
## Year_fct2018:WaterTemp                      0.019415 *  
## Year_fct2019:WaterTemp                      8.39e-05 ***
## StationCodeRD22:WaterTemp                   0.000832 ***
## StationCodeI80:WaterTemp                    0.829731    
## StationCodeLIS:WaterTemp                    6.85e-15 ***
## StationCodeSTTD:WaterTemp                   0.019208 *  
## StationCodeLIB:WaterTemp                    0.120121    
## StationCodeRVB:WaterTemp                    0.026229 *  
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.467471    
## StationCodeI80:Flow:FlowActionPeriodDuring  0.460348    
## StationCodeLIS:Flow:FlowActionPeriodDuring  4.30e-08 ***
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.103086    
## StationCodeLIB:Flow:FlowActionPeriodDuring  0.987816    
## StationCodeRVB:Flow:FlowActionPeriodDuring  0.778648    
## StationCodeRD22:Flow:FlowActionPeriodAfter  0.097598 .  
## StationCodeI80:Flow:FlowActionPeriodAfter   0.484364    
## StationCodeLIS:Flow:FlowActionPeriodAfter   2.34e-08 ***
## StationCodeSTTD:Flow:FlowActionPeriodAfter  1.22e-09 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter   1.50e-06 ***
## StationCodeRVB:Flow:FlowActionPeriodAfter   7.55e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(DOY) 3.848  3.988 38.09  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.885   Deviance explained = 88.8%
## -REML = 1434.9  Scale est. = 0.11336   n = 3478
appraise(m_gam_q_fa_wt)

k.check(m_gam_q_fa_wt)
##        k'     edf   k-index p-value
## s(DOY)  4 3.84787 0.9605145   0.015
draw(m_gam_q_fa_wt, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gam_q_fa_wt, pages = 1, all.terms = TRUE)

acf(residuals(m_gam_q_fa_wt))

Box.test(residuals(m_gam_q_fa_wt), lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  residuals(m_gam_q_fa_wt)
## X-squared = 6772.7, df = 20, p-value < 2.2e-16

Model with AR() correlation structure

Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.

# Define model formula as an object
f_gam_q_fa_wt <- as.formula(
  "Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)"
)

# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_fa_wt <- gamm(
  f_gam_q_fa_wt, 
  data = df_chla_wt_q,
  method = "REML"
)

m_gamm_q_fa_wt_ar1 <- gamm(
  f_gam_q_fa_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_fa_wt_ar2 <- gamm(
  f_gam_q_fa_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
  method = "REML"
)

m_gamm_q_fa_wt_ar3 <- gamm(
  f_gam_q_fa_wt, 
  data = df_chla_wt_q, 
  correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
  method = "REML"
)

# Compare models
anova(
  m_gamm_q_fa_wt$lme, 
  m_gamm_q_fa_wt_ar1$lme, 
  m_gamm_q_fa_wt_ar2$lme, 
  m_gamm_q_fa_wt_ar3$lme
)
##                        Model df       AIC       BIC    logLik   Test  L.Ratio
## m_gamm_q_fa_wt$lme         1 84  3037.848  3552.798 -1434.924                
## m_gamm_q_fa_wt_ar1$lme     2 85 -1962.244 -1441.163  1066.122 1 vs 2 5002.092
## m_gamm_q_fa_wt_ar2$lme     3 86 -1985.590 -1458.380  1078.795 2 vs 3   25.347
## m_gamm_q_fa_wt_ar3$lme     4 87 -1983.981 -1450.641  1078.991 3 vs 4    0.391
##                        p-value
## m_gamm_q_fa_wt$lme            
## m_gamm_q_fa_wt_ar1$lme  <.0001
## m_gamm_q_fa_wt_ar2$lme  <.0001
## m_gamm_q_fa_wt_ar3$lme  0.5318

Again, it looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.

AR(2) Model

# Pull out the residuals and the GAM model
resid_q_fa_wt_ar2 <- residuals(m_gamm_q_fa_wt_ar2$lme, type = "normalized")
m_gamm_q_fa_wt_ar2_gam <- m_gamm_q_fa_wt_ar2$gam

summary(m_gamm_q_fa_wt_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + 
##     Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
## 
## Parametric coefficients:
##                                               Estimate Std. Error t value
## (Intercept)                                  8.796e+00  6.363e-01  13.824
## Year_fct2016                                 1.238e+00  6.561e-01   1.887
## Year_fct2017                                 1.791e-01  5.995e-01   0.299
## Year_fct2018                                -4.588e-01  6.028e-01  -0.761
## Year_fct2019                                -1.297e-01  5.750e-01  -0.226
## StationCodeRD22                             -1.322e+00  6.714e-01  -1.969
## StationCodeI80                               4.878e-01  6.551e-01   0.745
## StationCodeLIS                              -1.669e+00  6.586e-01  -2.534
## StationCodeSTTD                             -3.777e-01  6.810e-01  -0.555
## StationCodeLIB                              -3.107e-01  7.561e-01  -0.411
## StationCodeRVB                              -1.787e+00  8.259e-01  -2.164
## Flow                                        -2.257e-04  3.659e-04  -0.617
## FlowActionPeriodDuring                       7.161e-02  1.103e-01   0.649
## FlowActionPeriodAfter                        1.159e-01  1.356e-01   0.855
## WaterTemp                                    1.664e-02  2.240e-02   0.743
## Year_fct2016:StationCodeRD22                -6.248e-01  5.824e-01  -1.073
## Year_fct2017:StationCodeRD22                -6.087e-01  5.523e-01  -1.102
## Year_fct2018:StationCodeRD22                -3.091e-01  5.418e-01  -0.571
## Year_fct2019:StationCodeRD22                -1.262e-01  5.437e-01  -0.232
## Year_fct2016:StationCodeI80                  1.584e-01  5.929e-01   0.267
## Year_fct2017:StationCodeI80                 -2.583e-01  5.628e-01  -0.459
## Year_fct2018:StationCodeI80                 -5.353e-01  5.522e-01  -0.969
## Year_fct2019:StationCodeI80                  4.972e-01  5.542e-01   0.897
## Year_fct2016:StationCodeLIS                  6.262e-01  5.677e-01   1.103
## Year_fct2017:StationCodeLIS                  1.909e-01  5.499e-01   0.347
## Year_fct2018:StationCodeLIS                 -1.089e-01  5.368e-01  -0.203
## Year_fct2019:StationCodeLIS                  3.895e-01  5.428e-01   0.718
## Year_fct2016:StationCodeSTTD                 3.693e-01  5.864e-01   0.630
## Year_fct2017:StationCodeSTTD                -6.079e-01  5.920e-01  -1.027
## Year_fct2018:StationCodeSTTD                -9.043e-01  5.621e-01  -1.609
## Year_fct2019:StationCodeSTTD                -1.126e+00  5.536e-01  -2.034
## Year_fct2016:StationCodeLIB                  9.232e-01  5.719e-01   1.614
## Year_fct2017:StationCodeLIB                 -3.232e-01  5.559e-01  -0.581
## Year_fct2018:StationCodeLIB                 -2.101e+00  6.040e-01  -3.478
## Year_fct2019:StationCodeLIB                 -8.633e-01  6.117e-01  -1.411
## Year_fct2016:StationCodeRVB                  2.471e-01  5.819e-01   0.425
## Year_fct2017:StationCodeRVB                 -4.778e-01  5.461e-01  -0.875
## Year_fct2018:StationCodeRVB                 -3.624e-01  5.353e-01  -0.677
## Year_fct2019:StationCodeRVB                 -7.749e-01  5.446e-01  -1.423
## StationCodeRD22:Flow                         6.934e-04  1.181e-03   0.587
## StationCodeI80:Flow                          4.800e-03  1.591e-03   3.017
## StationCodeLIS:Flow                         -2.042e-04  9.037e-04  -0.226
## StationCodeSTTD:Flow                         1.321e-04  1.033e-03   0.128
## StationCodeLIB:Flow                          4.830e-04  3.730e-04   1.295
## StationCodeRVB:Flow                          2.179e-04  3.660e-04   0.595
## Flow:FlowActionPeriodDuring                  2.747e-05  3.769e-04   0.073
## Flow:FlowActionPeriodAfter                   6.981e-05  7.286e-04   0.096
## StationCodeRD22:FlowActionPeriodDuring       7.152e-02  1.515e-01   0.472
## StationCodeI80:FlowActionPeriodDuring        3.482e-01  1.549e-01   2.248
## StationCodeLIS:FlowActionPeriodDuring       -7.908e-02  1.334e-01  -0.593
## StationCodeSTTD:FlowActionPeriodDuring       9.413e-03  1.343e-01   0.070
## StationCodeLIB:FlowActionPeriodDuring       -4.343e-01  1.791e-01  -2.426
## StationCodeRVB:FlowActionPeriodDuring       -1.575e-01  1.846e-01  -0.854
## StationCodeRD22:FlowActionPeriodAfter        2.049e-01  2.013e-01   1.018
## StationCodeI80:FlowActionPeriodAfter         4.531e-01  2.040e-01   2.221
## StationCodeLIS:FlowActionPeriodAfter         1.002e-01  1.688e-01   0.594
## StationCodeSTTD:FlowActionPeriodAfter       -4.901e-01  1.705e-01  -2.874
## StationCodeLIB:FlowActionPeriodAfter        -4.527e-01  2.014e-01  -2.247
## StationCodeRVB:FlowActionPeriodAfter        -1.957e-01  2.062e-01  -0.949
## Year_fct2016:WaterTemp                      -5.317e-02  2.012e-02  -2.642
## Year_fct2017:WaterTemp                       4.584e-03  1.825e-02   0.251
## Year_fct2018:WaterTemp                       2.569e-02  1.940e-02   1.325
## Year_fct2019:WaterTemp                       8.629e-03  1.730e-02   0.499
## StationCodeRD22:WaterTemp                    6.962e-02  2.193e-02   3.175
## StationCodeI80:WaterTemp                    -8.944e-03  2.043e-02  -0.438
## StationCodeLIS:WaterTemp                     4.623e-02  2.165e-02   2.135
## StationCodeSTTD:WaterTemp                    5.668e-03  2.295e-02   0.247
## StationCodeLIB:WaterTemp                    -3.515e-02  2.818e-02  -1.248
## StationCodeRVB:WaterTemp                     1.683e-02  3.139e-02   0.536
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.121e-03  1.164e-03  -0.963
## StationCodeI80:Flow:FlowActionPeriodDuring  -5.969e-03  1.557e-03  -3.832
## StationCodeLIS:Flow:FlowActionPeriodDuring   4.993e-04  9.440e-04   0.529
## StationCodeSTTD:Flow:FlowActionPeriodDuring  7.451e-04  1.066e-03   0.699
## StationCodeLIB:Flow:FlowActionPeriodDuring  -2.590e-04  3.895e-04  -0.665
## StationCodeRVB:Flow:FlowActionPeriodDuring  -1.933e-05  3.772e-04  -0.051
## StationCodeRD22:Flow:FlowActionPeriodAfter  -9.655e-04  2.553e-03  -0.378
## StationCodeI80:Flow:FlowActionPeriodAfter   -7.980e-03  2.766e-03  -2.885
## StationCodeLIS:Flow:FlowActionPeriodAfter   -6.062e-04  1.275e-03  -0.475
## StationCodeSTTD:Flow:FlowActionPeriodAfter   4.649e-03  1.391e-03   3.343
## StationCodeLIB:Flow:FlowActionPeriodAfter   -1.939e-04  7.341e-04  -0.264
## StationCodeRVB:Flow:FlowActionPeriodAfter   -6.629e-05  7.287e-04  -0.091
##                                             Pr(>|t|)    
## (Intercept)                                  < 2e-16 ***
## Year_fct2016                                0.059184 .  
## Year_fct2017                                0.765131    
## Year_fct2018                                0.446654    
## Year_fct2019                                0.821571    
## StationCodeRD22                             0.049085 *  
## StationCodeI80                              0.456585    
## StationCodeLIS                              0.011330 *  
## StationCodeSTTD                             0.579142    
## StationCodeLIB                              0.681142    
## StationCodeRVB                              0.030529 *  
## Flow                                        0.537306    
## FlowActionPeriodDuring                      0.516207    
## FlowActionPeriodAfter                       0.392757    
## WaterTemp                                   0.457427    
## Year_fct2016:StationCodeRD22                0.283391    
## Year_fct2017:StationCodeRD22                0.270529    
## Year_fct2018:StationCodeRD22                0.568362    
## Year_fct2019:StationCodeRD22                0.816454    
## Year_fct2016:StationCodeI80                 0.789413    
## Year_fct2017:StationCodeI80                 0.646245    
## Year_fct2018:StationCodeI80                 0.332393    
## Year_fct2019:StationCodeI80                 0.369731    
## Year_fct2016:StationCodeLIS                 0.270033    
## Year_fct2017:StationCodeLIS                 0.728490    
## Year_fct2018:StationCodeLIS                 0.839244    
## Year_fct2019:StationCodeLIS                 0.473020    
## Year_fct2016:StationCodeSTTD                0.528925    
## Year_fct2017:StationCodeSTTD                0.304603    
## Year_fct2018:StationCodeSTTD                0.107788    
## Year_fct2019:StationCodeSTTD                0.042061 *  
## Year_fct2016:StationCodeLIB                 0.106552    
## Year_fct2017:StationCodeLIB                 0.560996    
## Year_fct2018:StationCodeLIB                 0.000511 ***
## Year_fct2019:StationCodeLIB                 0.158224    
## Year_fct2016:StationCodeRVB                 0.671099    
## Year_fct2017:StationCodeRVB                 0.381721    
## Year_fct2018:StationCodeRVB                 0.498377    
## Year_fct2019:StationCodeRVB                 0.154848    
## StationCodeRD22:Flow                        0.557060    
## StationCodeI80:Flow                         0.002576 ** 
## StationCodeLIS:Flow                         0.821271    
## StationCodeSTTD:Flow                        0.898262    
## StationCodeLIB:Flow                         0.195464    
## StationCodeRVB:Flow                         0.551666    
## Flow:FlowActionPeriodDuring                 0.941913    
## Flow:FlowActionPeriodAfter                  0.923674    
## StationCodeRD22:FlowActionPeriodDuring      0.636993    
## StationCodeI80:FlowActionPeriodDuring       0.024636 *  
## StationCodeLIS:FlowActionPeriodDuring       0.553410    
## StationCodeSTTD:FlowActionPeriodDuring      0.944114    
## StationCodeLIB:FlowActionPeriodDuring       0.015339 *  
## StationCodeRVB:FlowActionPeriodDuring       0.393441    
## StationCodeRD22:FlowActionPeriodAfter       0.308635    
## StationCodeI80:FlowActionPeriodAfter        0.026442 *  
## StationCodeLIS:FlowActionPeriodAfter        0.552818    
## StationCodeSTTD:FlowActionPeriodAfter       0.004072 ** 
## StationCodeLIB:FlowActionPeriodAfter        0.024681 *  
## StationCodeRVB:FlowActionPeriodAfter        0.342716    
## Year_fct2016:WaterTemp                      0.008277 ** 
## Year_fct2017:WaterTemp                      0.801701    
## Year_fct2018:WaterTemp                      0.185344    
## Year_fct2019:WaterTemp                      0.618037    
## StationCodeRD22:WaterTemp                   0.001511 ** 
## StationCodeI80:WaterTemp                    0.661513    
## StationCodeLIS:WaterTemp                    0.032853 *  
## StationCodeSTTD:WaterTemp                   0.804969    
## StationCodeLIB:WaterTemp                    0.212274    
## StationCodeRVB:WaterTemp                    0.591847    
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.335423    
## StationCodeI80:Flow:FlowActionPeriodDuring  0.000129 ***
## StationCodeLIS:Flow:FlowActionPeriodDuring  0.596880    
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.484804    
## StationCodeLIB:Flow:FlowActionPeriodDuring  0.506064    
## StationCodeRVB:Flow:FlowActionPeriodDuring  0.959130    
## StationCodeRD22:Flow:FlowActionPeriodAfter  0.705263    
## StationCodeI80:Flow:FlowActionPeriodAfter   0.003938 ** 
## StationCodeLIS:Flow:FlowActionPeriodAfter   0.634486    
## StationCodeSTTD:Flow:FlowActionPeriodAfter  0.000838 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter   0.791721    
## StationCodeRVB:Flow:FlowActionPeriodAfter   0.927530    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value    
## s(DOY) 3.623  3.623 7.084 2.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.797   
##   Scale est. = 0.27259   n = 3478
appraise(m_gamm_q_fa_wt_ar2_gam)

k.check(m_gamm_q_fa_wt_ar2_gam)
##        k'      edf   k-index p-value
## s(DOY)  4 3.622926 0.8021592       0
draw(m_gamm_q_fa_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)

plot(m_gamm_q_fa_wt_ar2_gam, pages = 1, all.terms = TRUE)

acf(resid_q_fa_wt_ar2)

Box.test(resid_q_fa_wt_ar2, lag = 20, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  resid_q_fa_wt_ar2
## X-squared = 30.263, df = 20, p-value = 0.0657

The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.

# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_fa_wt_ar2_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + 
##     Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
## 
## Parametric Terms:
##                                   df     F  p-value
## Year_fct                           4 2.098 0.078512
## StationCode                        6 4.328 0.000234
## Flow                               1 0.381 0.537306
## FlowActionPeriod                   2 0.366 0.693855
## WaterTemp                          1 0.552 0.457427
## Year_fct:StationCode              24 2.863 4.02e-06
## StationCode:Flow                   6 3.598 0.001470
## Flow:FlowActionPeriod              2 0.005 0.994855
## StationCode:FlowActionPeriod      12 5.902 2.95e-10
## Year_fct:WaterTemp                 4 4.121 0.002471
## StationCode:WaterTemp              6 5.271 2.04e-05
## StationCode:Flow:FlowActionPeriod 12 4.904 4.35e-08
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F  p-value
## s(DOY) 3.623  3.623 7.084 2.26e-05

Despite its complexity, we’ll use this model for the model selection process.

Model selection - GAM models

As a summary, here are the 6 GAM models we are comparing:

  • Model 1 - m_gamm_cat_ar2 - Categorical predictors only
    Formula: Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5)

  • Model 2 - m_gamm_cat_wt_ar2b - Categorical predictors adding water temperature
    Formula: Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * FlowActionPeriod + s(DOY, k = 5)

  • Model 3 - m_gamm_q_ar2b - Flow as continuous predictor
    Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, k = 5)

  • Model 4 - m_gamm_q_wt_ar2 - Flow and water temperature as continuous predictors
    Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)

  • Model 5 - m_gamm_q_fa_ar2 - Flow with flow action period
    Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5)

  • Model 6 - m_gamm_q_fa_wt_ar2 - Flow with flow action period and water temperature
    Formula: Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
# AIC values
df_gamm_aic <- 
  AIC(
    m_gamm_cat_ar2$lme,
    m_gamm_cat_wt_ar2b$lme,
    m_gamm_q_ar2b$lme,
    m_gamm_q_wt_ar2$lme,
    m_gamm_q_fa_ar2$lme,
    m_gamm_q_fa_wt_ar2$lme
  ) %>% 
  as_tibble(rownames = "gamm_model")

# BIC values
df_gamm_bic <- 
  BIC(
    m_gamm_cat_ar2$lme,
    m_gamm_cat_wt_ar2b$lme,
    m_gamm_q_ar2b$lme,
    m_gamm_q_wt_ar2$lme,
    m_gamm_q_fa_ar2$lme,
    m_gamm_q_fa_wt_ar2$lme
  ) %>% 
  as_tibble(rownames = "gamm_model")

# Combine AIC and BIC and display results
left_join(df_gamm_aic, df_gamm_bic, by = join_by(gamm_model, df))
## # A tibble: 6 × 4
##   gamm_model                df    AIC    BIC
##   <chr>                  <dbl>  <dbl>  <dbl>
## 1 m_gamm_cat_ar2$lme        62 -2189. -1809.
## 2 m_gamm_cat_wt_ar2b$lme    65 -2204. -1806.
## 3 m_gamm_q_ar2b$lme         47 -2198. -1909.
## 4 m_gamm_q_wt_ar2$lme       58 -2174. -1818.
## 5 m_gamm_q_fa_ar2$lme       75 -2006. -1546.
## 6 m_gamm_q_fa_wt_ar2$lme    86 -1986. -1458.

According to AIC, model 2 (categorical predictors adding water temperature) was the best model, while BIC preferred model 3 (flow as continuous predictor). We’ll look at these two models more closely.

m_gamm_cat_wt_ar2b %>% saveRDS(here("manuscript_synthesis/model/interim/chla_cat_wt_gamm_model.rds"))
m_gamm_q_ar2b %>% saveRDS(here("manuscript_synthesis/model/interim/chla_flow_gamm_model.rds"))
df_chla_wt_q %>% saveRDS(here("manuscript_synthesis/model/interim/chla_wt_flow_model_data.rds"))